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INTRODUCTION 


Three-dimensional flutter solutions in subsonic and supersonic flows are 
an essential part of automated structural design. The methods that have been 
developed for subsonic flow are based on the numerical solution of the singular 
integral equation relating the pressure loading (or any other related quantity) 
to the downwash induced by the motion of the surface. The commonly used method 
is to treat the pressure as a series of preselected loading functions with un- 
known coefficients which are to be determined by satisfying the boundary 
conditions. In supersonic flows, the most frequently used method is to solve 
for the velocity potential directly In terms of distributed sources over the 
surface, although the integral equation approach has also been used. There is a 
need for a new method for both subsonic and supersonic speeds which is more 
oriented towards inclusion in a computer system (such as is used in automated 
design) . 

The objectives of this work include the development of a method based on 
the aerodynamic element concept which i s to be used over the subsonic and super- 
sonic Mach number range and the determination of aerodynamic forces on oscil- 
lating surfaces which are compatible with the structural breakdown. 

There are three principal features in the current developmental work. 
Firstly, the concept of aerodynamic elements is adapted to unsteady aero- 
dynamics. Secondly, the downwash-velocity potential method is used (also 
referred to as the Integrated potential method) which results in simpler 
expressions than those resulting from the normally used downwash-pressure 
formulation. Lastly, the method developed is applicable to both subsonic and 
supersonic speeds. 

In this report, a numerical method is formulated for calculating the 
aerodynamic matrix and generalized airforce matrix for both planar and non- 
pi anar wings of arbitrary planforms in steady and unsteady subsonic and super- 
son i c f I ows . 

The method uses the concept of "aerodynamic elements" and the downwash- 
velocity potential relationship (in preference to the downwash-pressure 



relationship more commonly used by others). The downwash-velocity potential 
method, using this aerodynamic element concept, was applied to lifting surfaces 
at steady angles of attack in a Note (Ref. I) by Haviland. in this case, the 
method was applied to a very simple steady-state rectangular wing model in sub- 
sonic flow. It was shown that a wake strip could be constructed using the same 
elements and that it could be terminated a few chord lengths from the wing. 

The arrangement of surface elements was selected in an attempt to simulate the 
nested horse-shoe vortices used by Hedman (Ref. 2). Exact agreement with the 
latter method was demonstrated. The work reported here is an extension of this 
method applied to the oscillating case in subsonic and supersonic flows. The 
application of the present method to an oscillating rectangular wing in sub- 
sonic flows has already been reported in reference 3. 

The idea of using arbitrary surface elements in the steady case has been 
advanced by Woodward (ref. 4) and by Rubbert and Saaris (ref. 5). In these 
cases, researchers intended to represent curved surfaces, fillets, and arbitrary 
planforms in an optimum manner. The use of closed vortex loops surrounding 
surfaces of uniform potential by Rubbert and Saaris Is of special interest here 
since It is the steady-state limit of the velocity potential method described 
in this report. Unfortunately, the concept of vortex lines has not proved to 
be useful for unsteady compressible flows. However, by specifying velocity 
potential distributions over the surface elements (without wake effects) and 
adding velocity potential distributions In the wake separately as functions of 
the trailing edge values, the computations are simplified considerably over 
those resulting from methods using the downwash-pressure methods. Both of 
these approaches employ doublet distributions to formulate nonhomogeneous 
integral equations. Examples of the downwash-pressure method (also called the 
integrated pressure method) include the kernel function method of Watkins, 
et a I . (ref. 6), and the doublet- lattice method of Albano and Rodden (ref. 7). 
Other researchers, such as Jones (ref. 8), Stark (ref. 9), Houbolt (ref. 10), 
and Morino and Kuo (ref. II) In the subsonic case, Rodemich and Andrew (ref. 12) 
in the transonic case, and Donato and Huhn (ref. 13) In the supersonic case 
have used the velocity potential. Most formulations of the supersonic problem 
have used the direct source method, as opposed to the doublet formulation 
necessary in the subsonic case. However, the downwash-pressure methods have 
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been employed by Watkins et al. (ref. 14) using the kernel function approach, 
and an extension of the doub I et- lattice method has been suggested by Harder 
and Rodden Cref. 15). A comprehensive review of the state of art was pro- 
vided by Landahl and Stark (ref. 16). 

The method of solution proposed in this work is to select Npp downwash 
collocation points at which the complex normal flow velocity components are 
known in terms of amplitudes and phase of oscillatory motion, and then to 
select an equal number of velocity potential distributions (each multiplied 
by an unknown factor which is to be found). In the simplest app I i cat ion, 
there is a collocation point for each surface element and one for each strip 
of wake. There is one point correspond i ng to each velocity potential at the 
center of a surface element or on the trailing edge. 


In the proposed application, one would first specify the surface element 
geometry, the collocation points, and the unknowns in the velocity potential. 

I One would then select the aerodynamic elements which one feels to be most 
appropriate and supply the geometric data specified for these. if the computer 

1 

system provides sufficient flexibility, one would be able to have subroutines 
jin the library which could process the input data so as to solve a given 
type of problem repeatedly. The basic computer program would compute an 
jaerodynamic influence matrix using subsonic or supersonic aerodynamics (which- 
, ever Is applicable) and would give complex downwash at the collocation points 
I in terms of the unknowns in the velocity potential. In most cases, the in- 
verse of this matrix would be required which has the nature of an "aerodynamic 
stiffness matrix" and gives the solution for the unknown velocity potential 
in terms of the downwashes. The next step depends largely on the particular 
application. In structural applications, the forces at the structural nodes 
would be calculated, using subroutines which would supplement the aerodynamic 
element subroutines. In flutter applications, similar subroutines would be 
used to calculate the generalized forces. The main objective of this work is 
to develop an efficient method for calculating the aerodynamic influence matrix 
Since most available results from various other methods are given in the form 
of overall lift and moment derivatives (which can be readily calculated as 
generalized forc'es), these were calculated in the present work. 
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In order to evaluate the effect iveness of the method discussed here, various 
comparisons were made against results obtained by using other methods found In 
the literature. At present, the subsonic rectangular elements have been 
successfully applied not only to steady and to harmonically oscillating rectan- 
gular wings but also to swept wings with or without control surfaces and to 
T-tail configurations. The demonstration of supersonic rectangular elements has 
been restricted to planar, steady-state rectangular wings of various aspect 
ratios. This has also been compared with results of previous research. 

The mathematical formulation of the problem begins with the well-known 
three-dimensional linearized potential flow equation In a moving coordinate 
system. This equation is deriveable by a perturbation of the Eulerian momentum 
equations, the continuity equation, and the equation of state. An alternative 
way of deriving this equation Is by using the Galilean transormatlon of the 
acoustical equation. The following assumptions are made: 

1 . 1 nv I sc 1 d flow. 

2. Adiabatic flow. 

3. Irrotational flow (except for a certain prescribed region 
downstream of a body). 

4. No body forces. 

5. Equation of state for perfect gas. 

The boundary conditions on a thin aerodynamic surface can be stated as 
fol lows: 

1. Normal flow velocities at the collocation points (often referred 
to as the "downwash") are determined by the requirement that flow 
conforms to upper and lower surfaces. 

2. Normal flow velocities are equal on upper and lower surfaces, 

(I.e., no relative velocities or "breathing" of the two surfaces). 

3. Disturbances are felt after their cause. This eliminates parts of 
the solution and is sometimes derived from the condition that there 
is no radiation from infinity (Sommerfeld radiation condition). 
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The wake region streaming behind an aerodynamic surface is a sheet of 
vortices. It is treated as an additional boundary which can accomodate any 
normal flow and, therefore, any potential discontinuity but cannot sustain a 
pressure differential. Beyond this measure, no further steps are taken to 
assure that the Kutta condition is met (i.e., that the pressure goes to zero 
at the trailing edge). 

The basic equation defines an acoustical perturbation of a uniform flow 
and is appropriate for subsonic, transonic and supersonic flows. However, no 
further consideration has been given to the transonic case because it is not 
appropriate for flows over bodies of finite thickness- Rather, for such cases, 
it would be necessary to perturb the actual, nonuniform flow over such a body. 
Experience to date with the comparison of theoretical and experimental flutter 
results indicates that wings of practical thickness can be treated successfully 
by uniform flow perturbat ions except in the transonic speed range. 

Further discussion on the derivation of the equation, assumptions, and 
limitations involved can be found in references 6, 10, II, and 14. 

In the following three sections, the theoretical development is described, 
a method of solution is discussed applicable to rectangular elements, and 
computer results are presented. These cover subsonic rectangu. elements 
applied to rectangular and swept wings and to T-tail configurations. In 
addition, supersonic elements are applied to steady state rectangular wings. 

The next section presents a proposed extended method for subsonic and 
supersonic elements of polygonal plan form. Formally, the new method Is based 
on the development of contour integrals around the elements, as i s a I so used 
for the steady-state supersonic results. Previously, however, surface integrals 
had been used for the subsonic case. Similar contour integrals have been used 
by Jones fref. 17). It can be readily demonstrated that the new method would 
give identically the same results as the former for all in-plane calculations. 

Although the main effort reported here has been concentrated on obtaining 
results for comparison with vortex lattice or doublet methods, the need for 
irregular aerodynamic elements may dictate the development of improved velocity 
potential distributions in which finite vortices would be eliminated. One 
method for doing this was described by Mercer et al. (ref. 18). However, 
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there is some expectation that the polynomial expressions described in the 
section "Proposed Extended Method of Analysis" may prove adequate for this 
purpose. 

The final sections contain conclusions and recommendations. They are 
followed by three Appendices, covering the subsonic computer program, standard 
integrals, and a description of the modifications made to the computer program 
to enable it to handle the steady-state supersonic calculations. 

THEORETICAL DEVELOPMENT 
Statement of Problem 

The problem of the present investigation is to determine the aerodynamic 
forces on thin, flat plates undergoing harmonic oscillations in subsonic and 
supersonic flows. 


Coordinate System 

The basic coordinate system is x’, y, z, or jr with a uniform flow velocity 
V in the positive x-direction as shown In Figure I. Each surface has its own 
coordinate system x, s, n, with x parallel to x’, s as the other coordinate 
in the plane of the surface and n as the local normal. The latter makes an 
angle g with the z-axis, measured in a right-handed manner about x. 

The x’, y and z system locates the "receiving point" or "collocation 
point," sometimes given as the point k. The sending point jp is designated by 
the dummy coordinates n, C and by the local coordinates a, v, and is 
sometimes given as the point Z. The angle between the local normal v and the 
z-axis is y- Surrounding each point Jl is an element of surface area sometimes 
referred to as a "sending region." 


Governing Differential Equation 

For the propagation of small disturbances that must be satisfied by 
the velocity potentials, the linearized potential flow equation in the co- 
ordinate system shown in Figure I may be written in the form: 
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(I) 


£:/ + V + ^ + »!i 

3xW 3x^2 3y^ 3z^ 


where V and c are local unperturbed flow velocity and sound velocity, 
respectively. 

I(0t 

For a harmonically oscillating system with $ equivalent to <|)(fr)e , 

Eq. ([) becomes 



/ \ 3^(j) 

Vc2 / 3x»2 


a2(f> ^ 3^ 3ft 

3x' 2 3y2 932 


( 2 ) 


Throughout the remainder of the report, the oscillatory case will be 
considered only, and the symbols for the field quantities will refer to their 
complex amplitudes, so that, for example, the actual velocity potential is 

This also applies to pressure p, to density p, and to normal velocity 

w. 


Integral Equation 

Since it is very unlikely that an exact solution to Eq. (2) can be found, 
one approach to solving the problem Is the transformation of the governing 
differential equation and its associated boundary condition to an integral 
equation for which approximate solutions can be found. For compressible flow, 
the integral equation as derived on the basis of the velocity potential and 
given for example in reference 1 3 may be written as 

(Mir) = / A((.4p) dS (3) 

This integral equation relates the velocity potential differentia! across a 
thin aerodynamic surface or a wake, A4>, which is taken to be positive If It 
decreases in the positive v-directlon, to the velocity potential In the field. 

In Eq. (3), G(j|r,fp), represents Green’s function which appears here In 
derivative or ’’doublet” form. The normal flow velocity w(Jr), at the 
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collcx:ation point <r where the local normal is n can be obtained by different- 
iating Eq. (3) with respect to n which results in 


w^r) = = / A(|)(j(p)K (rxp)dS (4) 

n-»-0 S ^ 

where 

= (5) 

^ n,v-»"0 

The surface of integration, S, appearing in Eq. (3) includes the 
aerodynamic surfaces and their wakes. Elsewhere In the field, A<|) Is zero. 

Expression for Green’s function .- The Green’s function, G<jr,jp), appearing 
in Eqs. (3) and (5) can be written in a general form: 


where 


G(fr,p) = U(Tp)Gp(jrjfp) + U(T^)G^(}r,fp ) (6) 

-iwTp 

Gp^(jr, 5 p) = -e V4 ttR 

-|0)T 

G^Cir^p) = -e /4irR 

r 2 = _ gt)2 + g2(y _ ^)2 + g2(2 - ^)2 

= {R - M(x' - 5’)}/3^c 
= Time of transit of "retarded wave’ 

= {-R - M(x’ - 5’)/3^c 
= Time of transit of "advanced wave" 
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U(T„) 

U(x^) 

In the subsonic case, (M < I), Is negative and G consists of only one 
term, Gp(fr,^p). For supersonic flow, (M > I), both and are positive 
when R is real; therefore, disturbances are restricted to the region of the 
aft Mach cone from the disturbing point. (This is equivalent ot the statement 
made later that only disturbances in the forward Mach cone from a given 
receiving point are felt.) Also note that Xp and x^ are used in the same sense 
as in reference 14. In supersonic flow, the retarded wave arrives ahead of the 
advanced wave. 


=! 

■I 


! if Xp is real and positive 
0 otherwise 


I if is real and positive 


0 otherwise 


Expression for K, 

4 > 

The term defihed by Eq. (5) might be termed the "velocity potential 

kernel function." It can be simplified when the thin surface coordinate system 

is used because K. varies with the normal coordinates n and v only through the 
<P 

"acoustical distance" R. Thus noting that 


^ 

3v \ 8v / \ R 3R 


with a similar relationship existing for the derivative with respect to n, 

K. can be written as 
9 



After substitution of the Green’s function from Eq. (6) the following 
expressions result: 
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Tiftr^p) = - -5- IB. \ = cosCy - g) (8) 

8n V 3v / 

Ki^r^P) = -|^f <9) 

' 32 $2 an ' r2 

9Ki 

K2(irrip) - -B^R ^ (ID 


In the expression for T 2 , the term n(iip) denotes the n-coordlnate of the 
sending point jp, and v({r) denotes the v-coordinate of the receiving pointer. 

In actuality, both Tj and Ta are equivalent to terms given by Vivian and 
Andrew (ref. 19). However, the expressions for Ki and Ka are different from 
the corresponding expressions given by these authors, and. In particular, they 
are only singular when tr and fp coincide. 

The Green’s function, G, Introduced In Eg. (5), is often referred to as the 
”free-field Green’s function" or as the "fundamental solution." The term 
Green’s function Is retained here for brevity and convenience. 

The kernel function, K., introduced In Eq. (5), relates the unknown 

<P 

velocity potential differential across an aerodynamic surface or Its wake to 
the known normal flow velocity or "downwash." Because of this relationship, 
the method is called the "Downwash-VelocI ty Potential Method." Other authors, 
such as Allen and Sadler (ref. 20), have used the term "Integrated Potential 
Method . " 

The alternative methods of Watkins et a I Cref, 6) and of Albano and 
Rodden (ref. 7) lead to kernels which have more complex forms and less tract- 
able singularities. The assumption has been made that the present form results 
in numerical simplicity; this Is, however, a difficult point to substantiate. 
Certainly, the aerodynamic element concept does not preclude the use of such 
kernels which would relate the unknown pressure differential to the known 
normal flow velocity as Is noted in the following section. However, the use 
of global defining functions for the unknown pressures, as introduced by 
Watkins, et ai (ref. 6), would lead to unnecessary complication and would 
negate the Inherent simplicity of the aerodynamic element concept. 



METHOD OF SOLUTION 


» 


Introduction 

The integral equation shown In Eq. (3) is solved using an approximate 
method. The method assumes that the velocity potential distributions over 
selected "elements" on the aerodynamic surfaces are expressed in terms of 
unknown parameters. The velocity potential distribution of each element is 
determined by satisfying normal velocity boundary conditions at a set of col- 
location points located over the surface using the relationship given in 
Eq. (4). The computer program sets up simultaneous equations relating the 
velocity potentials to the known normal flow velocities at collocation points. 
Typically, there is one such point to each element. Solution of the simul- 
taneous equations and interpretation of the velocity potential parameters 
results in the required aerodynamic forces. Since there is a one-to-one 
correspondence between velocity potential and pressure distribution, any system 
of calculation based on pressure distributiort can be performed In terms of 
velocity potentials. 


Integration by Discrete Elements 


In the downwash-i ntegra I equation as expressed in Eq. (4), the velocity 
potential distribution Acp^p) is the only unknown function within the 
equation. The integration of Eq. (4) is required for each collocation point 
k = I, . . ., Npp and can be performed by integrating separately over the 
appropriate sending regions relating to the unknown velocity potential 
amplitude coefficients (Z = I, . . ., N^p) . Thus, a set of simultaneous 
equations is formed which can be written as 


w 


k 


V 


'^RP 

Z=l ^ 


(12) 


where the matrix ^ is the "aerodynamic influence matrix." The principle 
concern of this report is the determination of the matrix A, . Once this 

K ^ 

matrix is known, it is a fairly stra ightforward procedure to obtain the 
perturbation velocity potential differentials A^’s for all the collocation 
points and to calculate the generalized airforce matrix. 
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Aerodynamic Elements 


In calculating the aerodynamic Influence matrix, the concept of "aerodynamic 
elements" Is employed. This means that the user Is given a catalog of sub- 
programs capable of handling elements of different shape and different Mach 
number range. For an experimental program such as reported In this work, this 
concept also provides a convenient way of evaluating various mathematical 
formulations of the same problem. The concept Is very similar to that used In 
finite element method In structural analysis In which there Is a basic program 
containing sufficient subprograms to handle a wide variety of structural 
elements which may be required. 

Basic form .- Expressions for the terms of the aerodynamic Influence 
coefficient ^ given In Eq. (12) can be obtained by breaking the surfaces up 
Into elements (the unknown velocity potential differential being associated 
with the element A). To facilitate the handling of wake effects and to pro- 
vide for cases where the velocity potential differential is not constant over 
an element. It Is convenient to define a corresponding point in each wake 
element A*. Thus, the velocity potential differential can be written in the 
form: 

At(A*,C,a) = At^ f Ol*,C,<T)exp{-iwCC(ii*) - C(Jl)II/V} (13) 

In the above expression, f Is a shape factor. However, even If this Is 
uniform, the exponential term assures that the potential at the point Jl* 
conforms to the expression relating the velocity potential In the wake to the 
velocity potential near trailing edge. In order to standardize the terminology 
for on-surface and wake elements, the notations £ and I* are used throughout, 
but it Is to be understood that they become Identical for on-surface elements. 

The system of points and regions is illustrated in Fig. 2, When Eq. (13) Is 
substituted Into Eq. (4), and the result is redefined as In Eq. (12), it Is 
found that the aerodynamic Influence matrix can be expressed in the form: 

A^ ^ = Z(£*)k^(k,£,£*) (14) 

which reduces for on-surface elements to 

\,£ " ( 15 ) 
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elements on wing 
designated by (£) 


Fig. 2. Designation of elements 
for Wing and Wake 
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where 


k^(k,£,£*) = // K^(r,5P)fCJt*,5,cj) exp{- C^a*) - CU)D}d^da (16) 

tn Eq. (16), the area Integration is restricted to the element JL*. 


Subsonic Aerodynamic elements .- Having performed the indicated 
differentiations in Eqs. (9) and (1.1), the expressions for Ki(jr,fp) and K 2 (ifr,{p) 
for subsonic flows are as follows: 


Ki (|r,p) 


U(T,) 


( 


Rc / 


Gp(ir,jp) 


CI7) 


K2dr,jp) 


U(T,) 


/ 33*^ 5ig)3^ 

\ R2 Rc 


cz) 


(18) 


where U^^p) and Gp([r,|p) are shown in Eq. (6). Using the expression in Eq. (7) 
and substituting the resulting expression for K^^r,|p) Into Eq. (16), k^(k,A,£^) 
becomes 


where 


k^(k,Jt,A*) = {cos (y - g)ki(k,A*) + k 2 (k,Jl,Jl*) 

io)L(k,£,£*) 

V 


exp ) - 


ki(k,£*) = - 


4tt 


/ / 


^ + Iw L 


R3(k) 


R2(k) 


exp 


lU) 


B^c 


CR(k) - R(k,£*) + M5 - M (£*)1 


X f(£*,5,a) dC da 


k2(k,£,£.^) =::^ // + 


£* 


R5(k) cR^(k) c2 R3(k) 


exp 


loi 


CR(k) - R(k,£*) + - M?(£*)I! 


X n(C,a)f (£*,C,a) d^ da 


(19) 


(20) 


( 21 ) 
15 .. 



(22) 


m2. m 

L(k,A,il*) = 5CA*) - + — UU*) - 5Ck)} + — R(k,jl*) 

3 ^ 32 

In the above expressions, R(k,Jl) Is R derived for the points k and A, 
while R(k) is R derived for the points k andfp. Also, n(?,a) is the n-coordinate 
of the point C,a on the sending element. 

Mid point constant potential (MPCP) rectangular elements ; The simplest 
derivation of the integrals in Eq. (19) can be obtained by insuring that the 
point A* is essentially the mid-point of the element A* and then evaluating the 
integrand at that point. Assuming a uniform potential distribution with 
f(A*,C,cr) equal to unity. 


ki(k,A*) 


A(A*) > 3^ ^ \j^ I 

4ir (R3(k,A*) c R2(k,jl*) 


(23) 


k2(k,A*) 


-A(A*)v(k)n(A*) 

4ir 


33 ^ 3i4o3^ 

R5(k,A*) cR‘t(k,A*) c^R3(k,A*) 

where A(A*) is the area of the element A*. For the remaining subsonic aero- 
dynamic elements, the same out-of -plane term as shown in Eq. (24) will be used. 

Zero order constant potential (ZOCP) rectangular elements ; In reference 
I, the same expressions were used as for the MPCP rectangular elements In the 
first evaluation of the In-plane steady state case. It was found that they did 
not provide correct results. Therefore, k^(k,A*) was rederived for a rectangular 
element. Since only the steady-state term, independent of the radial frequency 
0 ), was involved, it Is referred to as the zero order element. The resulting 
expression for kj(k,A*) can be written as: 

kj(k,A*) = kee(k,A*) (25) 

4tt cR^(k,A*) 

where 


(24) 
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d? da 


k = ^ / / J 

^ rectangle Jl* R^Ck) 





(26) 


Fs 2 (k,C»cr) can be obtained using the standard integral formula listed in 
Appendix B. 


SS 4Trv(k) 


tan 


-1 v(k)R(k) 
X(k)S(k) 


R(k) =7x^(k) + 3^S^(k) + 3^v^(k) 
X(k) = C - ?(k) 


C27) 

(28) 

(29) 


S( k) = a - a( k) (30) 

The terms 5^p^(5,*) etc. define the coordinates of the edges of the 
f"®ctangle making up the region and are shown in Fig. 3. When the receiving 
point k is coplanar with the sending region £*, v(k) is zero, and one must be 
careful to program the limiting equations accordingly. 

First order constant potential (FOCP) rectangular elements ; In order to 
improve the accuracy of the results, the complete integrand of Eq. (20) was 
expanded to the first order term in to, and then integrated. The result being: 

ki(k,)l*) = k25(k,)l*) I I R(k,i!,*) + M^Ol*) - MC(k) 

where k^^ is shown in Eq. (26) and 


]! 


+ — k^J3Ck,ll*) 


(31) 


k^stk,**) 


= !_/ ; 

4lT 

rectangle ji* 


{X(k)/R3(k)} d^ da 



^FOR^^*^ ) ^IN 


(32) 
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Fig. 3. Constant Potential Elements 




(33) 


logg{R(k) + BS(k>} 

Numerical difficulties will occur whenever a collocation point k is near 
a corner of the rectangular element I*. This situation should thus be avoided. 

First order zero pressure (FOZP) rectangular elements; If the potential 
Is varied in such a way that the distributed pressure remains zero, as shown 
in Fig. 4, the stepped potential of Fig. 3 is avoided, and the derivation 
of the forces acting on the elements is simplified. The modification consists 
of using the shape factor 

f(A*,C,CT) = exp {--iwCe “ (34) 


In Eq. (20) and the resulting expression for ki(k,Jl*) is written as: 


ki(k,il*) = jl + — [MR(k,£*) + - C(k)D 

i 


+ itok^g(k,£*)/V 
(35) 


This Is slightly different in form from that given in Eq. (31). 


Supersonic aerodynamic elements .- As a first step in the development of 
supersonic rectangular elements, the in-plane case is considered. Much of the 
subsonic formulation is directly applicable to the supersonic formulations. 
Completing the differential indicated in Eq. (9), the expression for Ki(tr,ip) 
in supersonic flow becomes 






1 


3 1 

c: 

H 

73 

lb) 


4ir 

-R 

3R 1 

I R ) 

c 

R2 J 

1 


3 1 

1 1 

1 ib) 

u(T^n 

47T 

_R 

3R 1 

R ' 

1 c 

R2 J 


I h) 

3^c 


<0) 


{R-M(x’-C)} 


{-R-M(x’-U> 


(36) 


Now Eq. (16) can be written in the form 
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Fig. 4. Zero Pressure Elements 


I 





= {kjp^(k,Jl*) + k,^(k,A*)} exp{- LCk,Jl,A*)} 


lA 


( 37 ) 


where 




4ir 


element A*[_R 9 R ( R 


exp 


- 

_ B^c 


{R(k) + MC - M^CA*)} fa*,K,o) da 


) ) iw U(Tp)1 

] 


(38) 


=k .! ! 




j^<vj 

id) U(t.) 
+ ^ 


8 R 

1 R 1 

c R^ 


id) 




X exp I- {-R(k) + - M?(£*) f(Jt*,C,a) d? da 

B^c J 


m2 

L(k,A,£*) = C(Jl*) - C(A) + " 5 (k)} 

B^ 

For a uniform potential distribution with f(A^,C,a) equal to unity. 


(39) 


(40) 


k,„(k,A*) = J kcc(k,A*)i:i + — {R(k,Jl*) + M^U*) - MC(k)}] 

IK j bb 


kus'k,**’! expj- 


i(i)R( k,£*) 
B^c 


( 41 ) 


J 


k,.(k,a*) = ^ k<,^(k,8.*) Cl + — {-R(k,il*) + MCCA*) - Me(k)>l 

. iwM , ,, ) i !a»RCk,)l*) ) 

+ -^r- k,,o(k,)l*) } exp{ = > 


c 'US 


6^c 


( 42 ) 


Thus, the expression for k.(k,J^-*) becomes 
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k^CI<,£*) = 21<3g (k,£^^) 


I + — r?CA*) - 5(k: 
e^c L 


X cos 


^ sin 


32c 




32c 


( 43 ) 


where 


k (k £*) = f f — ^ ^ 

^ 4tt J J R(k) 9R(k) 

element A* 


R(k) 


d? da (44) 


k^^(k, 


-' a 

element 


^ dc da (45) 


In-plane steady supersonic (IPSS) rectangular elements ; In the steady-state 
case, Eq. (43) reduces to 


k^(k,Jl*) = 2 kec (k,jl*) (46) 

(p OO 

where k^^ (k,Jl*) is the expression given in Eq, (44). If one makes the follow- 
ing substitutions: 


U = U(Tp,T^) 

X = X(k) = c - C(k) 

S = S(k) = 0 - a(k) 

(47) 

N = N(k) = V - v(k) 

a2 = f^2 - I 

R2 = R2(k) = X^ - a2s2 - 02 n 2 
then, in the planar case, where v(k) -»■ 0, 
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(48) 


ai ^ / y \ ^ /y. \ 

R 3R\Rr 3N2'VR/ 


but 

(_a2il.+ lL + li_W>J\ = 0 
' ax^ as^ aN^/'R/ 


Therefore, 


where 




4‘rr 


// 

e J ement 


4ir 


■ // 

element 




Using Green's theorem (ref. 21) 

kgg(k,Jl*) = dX + Q dS] 



The cutoff function In the above equation may be written 


I dx: 

In the form: 


U = U(-a|s| - X) 


with 


3U _ I 9U 

ax ~ a as 


(49) 


(50) 


(51 ) 


(52) 


(53) 
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Derivatives of the step function give rise to the Impulse function, otherwise 
called the Dirac delta function. For example. 


lx U(x) = «(x) =0 (X f 0) 

lx UCx) = 6(x) = “ Cx = 0) 

A useful integral property of this delta function Is 

/ f(x)6(x) dx = f(0) 


(54) 


(55) 


The next step Is to make use of these properties and to perform the line 
Integration indicated In Eq. (51) for complete and partial rectangular elements 
as sketched In Fig. 5. a. The figure shows an example of a rectangular wing 
broken up In 3 x 3 rectangular elements. Mach lines from one receiving point 
are shown to illustrate the complete and partial elements inside the aft Mach 
lines. Thus, the system of points and regions for steady supersonic calcula- 
tions Is slightly different from the subsonic case as shown In Fig. 2. The 
collocation points lie on the forward edge of each element. There Is no need 
to consider the wake element in this case. The particular arrangement of 
collocation points used in these calculations Is shown In Fig. 5. a. it is 
different from that used in the subsonic case but was selected for convenience, 
and does not, by any means, represent the final method recommended. 


The integration along the periphery of a complete element (see Fig. 5.b) 

gives 


kgg(k,A*) 


4tt 



R 


R 

+ 

R 


B 

XS 

A 

XS 

c 

XS 

D 


(56) 


For partial elements such as No. 2 in Fig. 5. a which are cut by Mach lines, 
the path of Integration is taken as indicated in Fig. 5.c. In this case. 


kss(k,i^*) 


-a2 , B/ I 8U X V ., 

4tt \ aR as r3 ) ^ 


+ 


47T 




( 57 ) 
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For the elements No. 3 and 4 in Fig. 5. a, the line Integrations similar to 
those used for elements No. I and No. 2 are no longer applicable because of 
the singularity problem when XCk} = 0. 

However, the nature of this singularity Is well known because It results 
from the two-dimensional flow solution relating pressure at a point to local 
angle of attack at that point. According to this, the two-dimensional com- 
ponent of normal velocity under a lift line is character Ized by a Dirac delta 
function . 

Therefore, the solution for the cross-section of a two-dimensional 
supersonic wing consists of a stepped section. However, the section Is known 
to be character Ized by a straight line at a uniform angle of attack. Either 
one must abandon the uniform potential element, or one must use some suitable 
smoothing process in which a smooth curve Is found which is a best fit to the 
stepped profile. The latter approach Is used here. 

The aerodynamic matrix is obtained in two parts. The first part gives the 
three-dimensional relationship along the lines described already. The second 
gives the two-dimensional component and in order to obtain this part correctly, 
the arrangement of collocation points shown in Figure 5 was selected. It will 
be noted that each collocation point is on the forward edge of a corresponding 
element (also a lift line because this Is where the velocity potential jumps). 
Thus, the two-dimensional contribution to the aerodynamic matrix can be ob- 
tained by Inverting the two-dimensional expression for lift due to angle of 
attack. This leads to the following procedure for the contributions of elements 
No. 3 and 4 to the downwash at collocation point No. 4. 

The lift coefficient per unit span of a two-dimensional lifting flat 
p I ate i s 


C = L = 4g 

^ pV^Ax /m^ - i" 

Lift is related to velocity potential as follows: 

L ^ 

pV^Ax "J^^x 


C58) 


(59) 
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Downwash is related to velocity potential by 


A(|) /-i/f/|2 - I \ 

V” ^ ^ y - a 

Comparing this with Eq. (12), it may be seen that 


C60) 


A 


k,k 


- 

2Ax 


= -K 


(61) 


Thus, the contribution to downwash at collocation point No. 4 (= k) from 
element No. 4, as in Fig. 5.d, is -K, while the contribution from element No. 3 

consists of the sum of the values K (because collocation No. 4 Is on the near 

R 1 

edge of element No. 3) and I . This scheme represents a process of 

i E.F 

averaging out the downwash over the element resulting from the physical 
consideration that the downwash Is continuously distributed over the surface 
rather than concentrated along a line. 

It will be noted that the two-dimensional part of the aerodynamic matrix 
has -K on the diagonal, with +K at certain other elements. 

The listing of the additional subroutines used in steady supersonic 
calculations and modifications made In subroutine DMATR are included in 
Appendix C. 


Calculation of Aerodynamic Forces 

Basic approach .- Let the downwash w^ be the sum of the contributions from 
all of the complex amplitudes q. of generalized coordinate. Then the normal 
flow velocity at each collocation point, W|^, can be expressed in the complex 
form: 


w 


k 


V 


N 

I 

i = l 



N 

I 

i = l 






(62) 


where h. . and a; . are the amplitudes of displacement and local angle of 
I , K ^ ^ 

attack, respectively, at each col location point k for unit amplitude of the 
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generalized coordinate q., and N is the number of modes. 
Also, let 


N 

= ,1, ''i 


then 



'RP 

I 

k=l 



(63) 


(64) 


Thus, given the inverse of the aerodynamic Influence matrix A. and the modal 

^ p iC 

components of downwash W. . , the modal components of the velocity potential can 

i , K 

be calculated. 


Generalized forces .- The generalized forces Q. can be expressed in the 

form: 

N 

1 }! d( i ,x,s)AP( j,x,s)qj dx ds (65) 

j=l 


where d(i,x,s) and AP(j,x,s) are the modal deflection and net pressure at 
point (x,s), respectively. Eq. (65) can be put in the form: 


N 

Q = -pV2 I 
1 = 1 


C. . q . 

I »J J 


( 66 ) 


where p is the density of air. 

The flutter airforce matrix, C. j, is given as 


C- • = -I // d(i,x,s) dx ds 

'»J pV 2 


(67) 


Note that A(J) and AP are given in the same sense and are thus positive when the 
algebraic value is highest on the negative (-n) side of the airfoil. A positive 
AP then results in a force in the plus n direction. It can be shown that the 
pressure differential and velocity potential differential are related by the 
fol lowing equation : 
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AP(J,x,s) _ 
pV^ 

where f(Jl,x,s) is the shape factor factor mentioned in Eq. (13). Eq. (67) now 
can be expressed in terms of velocity potentials by the above relationship. 

On integration of Eq. (68), 


( 68 ) 


A(|>f (x,s) 
V 


-exp 



AP(A,s) 

pV^ 


exp 



(69) 


The velocity potential in the wake is related to the velocity potential near the 
trailing edge by the requirement that the pressure differential sustained by the 
wake Is zero. Thus, if A(f^^ is the velocity potential differential at x-p^ near 
the trailing edge, then 

A<J)f (x,s) = exp I - (x - Xj^) I (70) 

This can be derived readily from Eq. (69). 

Flutter airforce matrix for constant potential case .- in this case, the 

velocity potential varies as shown In Fig. 3, leading to lift lines at the 

front and rear of the rectangular element. The forward lift is negative (i.e., 

A(f): p 

in the negative n-d I rect ion ) . When the velocity potentials are solved 

from Eq. (64), the rectangular element moves as a rigid body so that 


d( i ,x,s) 




(71 ) 


The contribution to the airforce matrix is obtained from Eqs. (67) and (68) 
as fol lows: 



II [h 

element Z 


+ {x - x(£)}a 




9f (£,x,s) 
9x 


Id) 

V 


f (A,X,S) 


J dx ds I 


A(f 




(72) 
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Under the assumption that the velocity potential Is constant over the rectangle 
and that the point is at its center, the contribution to C. . for the surface 
e I ements becomes * 


AC. . -A(Jl) 

I * j 





(73) 


The velocity potential on each element in the wake is given by Eq. (70). 
The strength of the lifting line at the leading edge of the wake is thus equal 
to the velocity potential in the first element as follows: 


,s} = A(|)(£*) exp 


JtO 




- x(£*) + ^ Ax(£*)H 


(74) 


and the contribution to C. . is 

1 .J 


AC.^. =As(£*) + [:x,q,(£-) -x(.-):a.^^} 

(75) 

Att n* ( tW , ) 

^ — exp |- v~ 1 Ax(A*)H | 


where Xp^^Oi,*) is the x-coordinate of the leading edge of the wake strip, and 
Ax(£*) Is the chord of the first rectangular element in the wake, as shown 
In Fig. 3. If the wake is terminated at a finite distance (as in Fig. 3), 
there is a lift line in the wake which has negligible effect and is therefore 
ignored . 

Flutter airforce matrix for zero pressure case .- The proper interpretation 
of the pressures in the wake is open to question when constant potential 
elements are used. No such questions arise, however, when the zero pressure 
element is used. Here the velocity potential is assumed to vary across a 
rectangular element in such a way that the pressure is zero except under lift 
lines at forward and rearward edges. Thus, in a region influenced by the 
velocity potential A<|>(Jl), the shape factor Is expressed as 
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f (£,x,s) = exp ^ Cx - xCJDj | 


(76) 


The contribution to the flutter airforce matrix in this case is obtained by 
applying Eq. (72) and using the shape factor f(^-,xSs) varying as in Eq, (76). 


= Asia) + {XpQp(£) - xU)} a.^^H 

X exp j^- ~ {XpQp(^) - x(A)}J 

’ “i,A^ 

X exp ^ I — 


(77) 


The lift at the leading edge of the wake strip is similar to the lift at 
the leading edge of a rectangular region in Eq. (77), so that 



= ASCII*) Ch,^^, + {Xpo^(£«) - xOl*)} 
X exp {XrOr''^*' ■ 


(78) 


COMPUTED RESULTS AND DISCUSSIONS 


Computer Program 

A computer program was written in FORTRAN IV language to perform the 
calculations outlined in the section on METHOD OF SOLUTION. The description 
of the program, input instructions, and sample problems along with the list- 
ing of the program are included In Appendix A. 

To date, four versions of aerodynamic elements all using rectangular 
elements have been formulated. They are summarized as follows: 
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Subroutines 


Type of Elements 


AC0N5, AC0N6 and 
CC0N5, CC0N6 

AC0N7, AC0N8 and 
CC0N5, CC0N7 

ACONI I , AC0NI2 and 
CCONII, CC0NI2 

AC0N3, AC0N4 and 
CCON I 1 , CCON 1 2 


Zero Order Constant Potential (ZOCP) 
Rectangular Elements 

First Order Constant Potential (FOCP) 
Rectangular Elements 

First Order Zero Pressure (FOZP) 
Rectangular Elements 

In-Plane Steady Supersonic (IPSS) 
Rectangular Elements 


Scope of Computations 

A total of eight planforms were considered and are shown in Figures 6 to 
13 with the details of their geometry. Table I gives the details of the cases 
for which calculations were made. Subsonic calculations including cases I 
through 3 are presented In reference 3. 

Discreti zation of Surfaces 

The breaking up or division of a surface into rectangular elements was 
accomplished by defining a chordwise strip of wing, locating its centerline, 
dividing the length into a specified number of chordwise elements, and laying 
these out by shifting one quarter of the chord of an element back from the 
leading edge. Collocation points were then located at their centers. These 
were extended into the wake with approximate! y equal chordwise dimensions and 
terminated about five to ten wing chord lengths back. This method of locating 
elements is illustrated in Fig. 14 for a rectangular wing and is compared with 
a similar wing divided up into winglets with a lift line at the quarter chord 
of each and a collocation point at the three-quarter chord as employed by 
Hedman (ref. 2). 

Only the half wing is shown divided up in Fig. 14, because the symmetry 
logic of the computer program takes care of the effects of the other half. 
Other planforms used in the present calculations also have one or two planes 
of symmetry and these symmetries were used whenever applicable in an attempt 
to save computational effort. 
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Fig. 



Rectangular Wing A.R. = 2 (Planform No. 1) 












Mach No . 

1.2 
1.4 
1 . 6 
1.8 


j Mach No. = 1.3 

» Aspect Ratio = 1.2039 

M— * 

Mach No. = 1.3 
Aspect Ratio = 2.0 


Mach No. = 1.3 
Aspect Ratio = 4.0 


Fig. 13 Rectangular Wings for Supersonic Calculations 
(Planform No. 8) 
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Table I 


Calculations Performed 


Case 

No. 

Planform 
Number j 

Mach 

Number 

Reduced 

Frequency 

Derivatives 

Remarks 

■ 

1 . Rectangu 1 ar Wing 
A.R. = 2 

0.9 

0.5 


Evaluation of 
the effective 
wake length 

2 

1 . Rectangular Wing 
A.R. = 2 

0 

0.8, 0.9 

0, .2, .4, 
.6, .8, .9 

0, .15, .25, 
.35, .5 

^La'^Ma 

Comparison with 
Kernel function 
method and 
Doublet Lattice 
method 

3 

2. Rectangular Wing 
A.R. = 8 

0 

0 


Strip width 
variation, and 
Comparison with 
Kernel function 
method 


3. Tapered Wing, 
A.R. = 5, 

Taper Ratio = .5 
Sweep angle 
= 3.817 deg 

.15 

0 

^m 

^PM'^L 

^PM'^BM 
BMCP,PMCP 
See "Results 
and Discus- 
sion, Case 4 

Comparison with 
tabulation by 

Thomas and 
Wang (ref. 26) 

4. Untapered Wing 
A.R. = 5 
Sweep angle 
= 15.0 deg 

■ 

0 

5 

5 . Swept Wing with 
Partial span 
flap A.R. = 2, 
L.E. Sweep angle 
= 60 deg 

0, .7, 
0.8, .9 

0.5 

C. . 

1 .J 

(See 

Appendix A) 

Comparison with 
the results re- 
ported by Stark 
(ref. 24) 

6 

6. Rectangular 
T-ta i 1 

0, .25 

0.2 


Comparison with 
Clevenson & 
Leadbetter 
(ref. 28) and 
Kalman (ref. 23) 
et . a 1 . 

7 . Stark’s T-ta i 1 

0.8 

0.3, .45 

^m,n 

Comparison with 
Stark’s results 
(ref. 24) 
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Table I Continued 


Calculation Performed 


Case 

No. 

Planform 

Number 

Mach 

Number 

Reduced 

Frequency 

Derivatives 

Remarks 

7 

8. Rectangular Wings 
with various 
Aspect Ratios 

1.2, 1.4, 
1.6, 1.8 

0 

c, ,c». 

La' Ma 

(See "Results 
and Discus- . 
si on. Case 7) 

Comparison with 
Ref. 30 


8. Rectangular Wings 
with Various 
Aspect Ratios 

1 .3 

0 

c, ,c,. 

La' Ma 

Comparslon with 
Nel son, et al . 
(ref. 31) 
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Fig. 14. Comparison of Four^by^-^Four Arrays, 
Present Method and Hedman's fref,.2) 


Results and Discussions 


Case I; convergence of the resul ts w ith the wake length .- Case I 
consisted of a check into the convergence of the results when the number of 
chord lengths included in the wake is varied. This is expressed as the ratio 
F which Is equal to the length of the wake divided by the chord. Results 
shown in Table II were obtained for a four-by-four arrangement on a rectan- 
gular wing (planform No. I) using first order zero pressure rectangular 
elements. The reduced frequency was 0.5, and the Mach number was 0.9. 

(Reduced frequency Is defined as wb/V, where b is the half chord of the wing), 

TABLE 1 I 

Convergence of Calculated Results for Aspect Ratio 
Two Rectangular Wings, Oscillating about Mid-Chord at Reduced Frequency of 
0.5 and Mach Number 0.9, as Ratio Wake Length/Chord (=F) is Changed 


Ratio 

Lift 1 

Derivative 

Moment 

Derivative 

F 

Rea 1 

Imag i nary 

Rea 1 

1 mag i nary 

I 

3.5454 

0.5387 

.6146 

-1 .8636 

5 

3.5622 

0.5444 

.5986 

-1 .8542 

10 

3.5622 

0.5440 

.5992 

-1.8542 

20 

3.5623 

0.5435 

.5996 

- 1 . 8542 


In the remaining calculations, F was taken equal to ten (10) which Is 
well within one percent error according to Table II. 

Case 2: rectangular wing pitching about , midchord . - Case 2 was concerned 

with calculations on a rectangular wing (planform No. I) to compare with 
two different references. First, calculations were made to compare with the 
results of the kernel function method reported by Runyan and Woolston (ref. 22) 
for the zero Mach number case and for reduced frequencies ranging up to 0.9. 
Second/ calculations were made to cpmpare wrth the ^resu kts of the doublet lattice 
method reported by Kalman, Rodden, and Glesing (ref. 23) for Mach numbers of 
0.8 and 0.9, and for reduced frequencies up to 0.5. These results have been 
reported in reference 3. 
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Using the key to the symbols given in Fig. 15, the results of these 
comparisons can be seen In Figs. 15 to 18 for comparison with the kernel 
function method and In Figs. 19 to 22 for comparison with the doublet lattice 
method. In the first comparison, generally fair agreement is noted except 
that the real part of the moment derivative does not agree well at zero fre- 
quency. Such a disagreement was noted earlier (ref. I) when it was seen that 
the kernel function method and the downwash-ve locity potential method gave 
different centers of pressure on a high aspect ratio wing. Otherwise, agreement 
with the first order zero pressure elements is generally high. The constant 
potential elements have the next highest agreement and the zero order calcu- 
lations appear to have least agreement. At zero frequency, these three methods 
are of course identical . Improved agreement is also noted with the eight-by- 
eight arrangement as compared to the four-by-four arrangement. This indicates 
convergence as the number of collocation points (Npp) is increased. One result 
by Stark (ref. 24) is shown for comparison and Is seen to be in very close 
agreement with the results of Runyan and Woolston (ref, 22). In the second 
comparison, agreement at zero frequency is now good. This should be expected 
considering that the doublet lattice method gives the Hedman results at zero 
frequency. Agreement at a reduced frequency of 0.5 and a Mach number of 0.9 is 
good when first order zero pressure elements are used in an eight-by-elght array. 
Results given by Laschka (ref. 25) and one result by Stark (ref. 24) are also 
shown. 


Definitions of the derivatives shown in the figures are as follows: 

wc 

Reduced Frequency “ 


3L 


Lift Derivative “ 9^ 

Moment Derivative = ^ 4 ciSc 

where L is the lift in the sense of the negative n-axis, a is 
rotation about the mid chord, q is the dynamic pressure, S is 
M is the moment about the mid chord line, and c is the chord, 
are positive in the same sense. 


the angle, of 
the wing area. 
Both a and M 
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SOLID LINE- RUNYAN S WOOLSTON 

B-ZERO ORDER, 4x4 ARRAY 

A-FIRST ORDER, CONSTANT POTENTIAL, 4x4 ARRAY 

O-FIRST ORDER, ZERO PRESSURE, 4x4 ARRAY 

■ , A , 0-8x8 ARRAYS 



LIFT 

DERIVATIVE 

(REAL) 


1 ^ 


REDUCED FREQUENCY 

^0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 


Fig. 15. Lift Derivative (Heal Part) for Aspect Ratio Two Rectangular Wing Rotating 
about Mid**Chord at Zero Mach Nuober 



Fig. 16. Lift Derivative (Imaginary Part) for Aspect Ratio Two Rectangular Wing Rotating 
about Mid*-Chord at Zero Mach Number 



Fig. 17. Moment Derivative (Real Part) for Aspect Ratio Two Rectangular Wing Rotating 
about Mid-Chord at Zero Mach Number 





1.6 



SOLID LINE-KALMAN, RODDEN, ANDGIESING 
DOTTED LINE-LASCHK^. 

X-STARK 

OTHER SYMBOLS-SEE F IGURE 15 


1.0 F 


0 . 8 | 


MOMENT 

DERIVATIVE 

(REAL) 


0.6 f 


0.4 F 



REDUCED FREOUENCY 

OJ 02 ^"o!3 04 05 


Fig. 19 Moment Derivative (Real Part) for Aspect Ratio Two Rectangular 
Wing Rotating about Mid- Chord at Mach Number 0.8. 





Piy. 20 Moment Derivative (Imaginary Part) for Aspect Ratio 
Two Rectangular Viing Rotating about Mid-Chord at 
Mach Number 0.8. 








Pig. 22, Moment Derivative (Imaginary Part) for Aspect Ratio Two 

Rectangular Wing Rotating about Mid-Chord at Mach Number 0.9. 


Case 3: lift distribution of a rectangular wing; subsonic steady .- 

Sectional lift slope, aerodynamic center total lift, and moment derive 
tlves for an aspect ratio eight rectangular wing (planform No. 2) v/ere 
evaluated in this case. In order to examine the effect of the strip width 
variations, one set of runs was made with an evenly spaced strip width as 
shown in Fig. 23. Another set of runs was made with the arrangement having 
reduced strip width towards the tip as shown in Fig. 24. The lift distri- 
butions and aerodynamic centers are plotted in Fig. 25 and Fig. 26, respec- 
tively, with results of the kernel function method for comparison (obtained 
by private communication from E. C. Yates, NASA Langley Research Center), 

The sectional lift slope distribution in Fig. 25 shows that the present 
method gives lower values than the kernel function values, although the trend 
reverses towards the tip. Improved agreement was also noticed towards the 
tip when the arrangement of reduced strip width towards the tip, sketched in 
Fig. 24, was used. The aerodynamic centers calculated by the kernel function 
method using polynomial pressure modes are in closer agreement with the 
present results than those obtained by the kernel function method using the 
Glauert-Bi rnbaum series. This is shown In Fig. 25 with numerical data in 
Table 111. 

Case 4; tapered and swept wings; subsonic steady .- In this case, 
calculations were performed at the request of Langan and V/ang of NSRDC .<ref. 26) 
for the comparison of various numerical lifting surface theories with experi- 
mental results. Calculations were made for a tapered wing of aspect ratio five 
(planform No. 3) and for a swept wing of the same aspect ratio with no taper 
t planform No. 4). The arrangements of the aerodynamic elements used for these 
calculations are shown in Figs. 27 and 28. 

Figures 29 and 30, and Tables IV and V were supplied by Langan and 
Viang '(ref. 26). Tables IV and V present the wing loading characteristics for the 
tapered wing and the swept wing as computed by 15 programs and as obtained 
from the experimental results. The programs are listed In Tables IV and 
V and are Identified by the names of those who contributed the data. Mach 
number for which the calculations were made, lift coefficient Cj^, pitching 
moment center of pressure (PMCP) measured from the wing apex divided by the 
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Kernel Function Method - 40 Stations (4 ^ 4) , Glauert Bimbaum Series 

First Order, Zero Pressure, 8x8 Array, Evenly Spaced Strips 

First Order, Zero Pressure, 8x8 Array, Reduced Strip Width Towards the Tip 

Kernel Function Method - (6x6), Langi'.ey Kernel Function Program Using 
Polynomial Pressure Modes 



ROOT TIP 


Fig. 25 Comparison of Section Lift Slope For An Aspect Ratio 8 Rectangular Wing 
Mach No. a 0, k 0 


Kernel Function Method ~ 40 Stations (4 x 4) , Glaiiert-Birbaum Series 


First Order, Zero Pressure, 8x8 Array, Evenly Spaced Strips 

First Order, Zero Pressure, 8x8 Array, Reduced Strip Width Towards the Tip 



Fig. 26 Comparison of Aerodyncunic Center For An Aspect Ratio 8 Rectangular Wing 
Mach No. » 0, k = 0 


Table III 

Sectional Lift Slope and Aerodyncunic Centers, Total Lift and 
Moment for Aspect Ratio 8 Rectangular Wing with the variations 
in Station Spacing. Mach number =0, *= 0, constemt down- 

wash a » 1.0 


Aerodynamic 

Element 

Strip 

Spacing 

y 

” “ b 


Cm 



Cm 

First Order 
Zero Press. 

Even 

Spacing 

(See 

Fig. 23) 

.0625 

5.3187 

1.3194 

0.2481 

4.7226 

1.1474 

.1875 

5.2896 

1.3113 

0.2479 

.3125 

5.2253 

1.2931 

0.2475 

.4375 

5.1147 

1.2616 

0.2467 

.5625 

4.9359 

1.2102 

0.2452 

- 

.6875 

4.6426 

1.1248 

0.2423 

.8125 

4.1314 

0.9747 

0.2359 

.9375 

3.1128 

0.6943 

0.2198 

First Order 
Zero Press. 

^^djusted 

spacing 

with 

reduced 

strips 

towards 

4-1 r\ 

.1250 

5.3132 

1.3177 

0.2480 

4.6898 

1.1393 

.3750 

5,1936 

1.2844 

0.2473 

.5625 

4.9559 

1.2164 

0.2454 

.6875 

4.6588 

1.1300 

0.2426 

(see 

r>r "J A \ 

.78125 

4.2873 

1.0206 

0.2381 


.84375 

3.8688 

0.8976 

0.2320 


.90625 

3.2858 

0.7308 

0.2224 


.96875 

2.3371 

0.4809 

0.2058 
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Fig. 29 Spanwise Lift Distribution for 
From Langan and Wang^® 


Wing at 11.4® Angle of Attack 







o 

o 



Fig. 30 - Spanwise Distribution of Chordwise Location of Center of Pressure for the 

Tapered Wing at 11.4-Degree Angle of Attack From Langan and Wang26 










TABLE IV - Tapered Wing at 11.4-Degree Angle of Attack 



Hach 


PHCP 

8HCP 

Si 

NCM X NSM 

Tim# 

sec 

Computer 

Tullntus 

0. 

15 

0.817 

0.240 

0.425 

0.043 

420^*^ 

25 

CDC 6600 

Dulmovtts 



0.831 

0.240 

0.431 


10 X 15 

96 

IBM 360/75 

Hargaton-Ltmar 



0.831 

0.240 

0.432 


10 X 12 

142 

CDC 6700 

Gitsing 

0. 

15 

0.829 

0.237 


0.044 

11 X 17 

120C 

IBH 360/65 

Rubbtrt 

0.00 

0.798 

0.238 

0.420 

0.041 

10 X 15 

117 

CDC 6600 

Loptz A Sh«n 

0.00 

0.814 

0.238 


0.042 

8 X 17 

105C 

IBM 360/65 

HavlTaml 

0. 

15 

0.848 

0.242 

0.438 


8x8 

291 

CDC 6400 

Jordan 

0. 

00 

0.820 

0.211 


0.043 

2 X 15 

9 

UNIVAC 1108 

Lanwr 

0. 

00 

0.806 

0.239 

0.423 

0.041 

4 X 13 

98 

IBM 360/65 

Wldnall 

0. 

15 

0.812 

0.218 

0.431 


10 X 2 

63 

CDC 6700^^^ 

Bandler 

0. 

00 

0.807 

0.238 


0.042 

4x4 

B 

COC 6700^^^ 

Rowt 

0. 

15 

0.815 

0.239 

0.428 


4x6 

13 

CDC 6600 

Cunningham 

0.15 

0.839 

0.235 

0.424 


4 X 5 

14 

IBM 370/155 

Jacobs-Tsakonas 

0.00 

0.920 

0.264 

0.439 


5(b) 

40C 

COC 6600 

Lopez (KDchemann) 

0. 

15 

0.852 

0.242 


0.047 

16^*=) 

12C 

IBM 360/65 

Experiment 

0. 

15 

0.83 

0.28 

0.44 





Additional Computer Runs 

Closing 

0. 

00 

0.823 

0.240 


0.044 

11 X 17 • 

121C 

IBM 360/75 

Jordan 



0.808 

0.212 


0.042 

2 X 7 

5 

UNIVA 

1108 




0.819 

0.211 


0.043 

! . IS*") 

.10 






0.817 

0.213 


0.043 

2 X 31 

28 






0.811 

o.zn 


0.042 

3x15 

31 



Jordan 

0. 

00 

0.798 

0.213 


0.041 

4x15 

68 

UNIVAC 1108 

Cunningham 

0. 

15 

0.849 

0.233 

0.425 


3 X 3 

10 

IBM 370/155 

Lopez 

0. 

00 

0.845 

0.242 


0.046 


12C 

IBM 360/65 

Widnall 

0. 

15 

0.826 

0.217 

0.426 


10 X 3 

73 

CDC 6700 

Notes: 











(a) Tullnlus uses 420 horseshoe vortices with 6 chordwise and 6 spanwise functions relating 
their strengths. The system Is solved by a least squares method with 99 pivotal points. 

(b) Chordal modes. 

(c) Spanwise modes. 










(d) The kink at 

the wing centerline has been rounded off for these calculations. 



(e) Experimental results are examined In detail In 
Experimental Results." There Is a difference 
and those In the calculations. 

the section entitled. "Comparisons with 
between the wing tips In the experiment 

(f) Times for CDC 6700 are In terms of CDC 6400 time. 






From Langan and Wang 
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TABLE V - Uncambered Swept Wing 


Progrt* 

Hach 

‘^Lo 

PMCP * 

8MCP 

^01 

NCM X NSM 

Time 

sec 

Computer 

Tullnlus 

0.12 

3.891 

0.536 

0.450 

0.999 

420^*^ 

87 

CDC 6600 

Duimvits 

0.12 

3.987 

0.541 

0.455 


10 X 15 

104 

IBM 360/75 

Mirgason-Ltmir 

0.10 

4.008 

0.532 



10 X 12 

71 

COC 6700^*^ 

Glesing 

0.12 

3.973 

0.540 


1.053 

11 X 17 

123C 

IBM 360/65 

Rubbert^**^ 

0.00 

3.868 

0.535 

0.447 

0.966 

10 X 15 

35^(c) 

COC 6600 

Lopez & Shen 

0.00 

3.916 

0.536 


1.000 

8 X 17 

HOC 

IBM 360/65 

Havlltnd 

0.12 

4.021 

0.583 

0.463 


6 X 10 

463 

COC 6400 

Jordan 

0.00 

3.910 

0.531 


0.999 

2 X IS 

13 

UNIVAC 1108 

Lamar 

0.00 

3.874 

0.536 

0.447 

0.980 

8 X 11 

651 

IBM 360/65 

Uidnall 

0.12 

4.080 

0.502 

0.443 


10 X 2 

60 

CDC 6700^®^ 

Bandler 

0.00 

4.030 

0.528 


1.052 

4x4 

10 

CDC 6700^®^ 

Rowe 

0.12 

3.921 

0.534 

0.447 


8x9 

95 

CDC 6600 

Cunningham 

0.12 

4.091 

0.523 

0.444 


5x5 

19 

IBM 370/155 

Jacobs-Tsakonas 

0.00 

4.253 

0.553 

0.464 


5(«*) 

80C . 

CDC 6600 

Lopez 

0.12 

4.065 

0.550 


1.109 

16 <*J 





Experiment | 

Pressure 

0.12 

3.89 


0.462 





Force 

0.12 

3.85 


0.466 





Additional Computer Runs 

Marga son- Lamar 

0.00 

3.996 

0.532 

0.457 


10 X 12 

144 

CDC 6700^®) 

Glesing 

0.00 

3.956 

0.540 


1.044 

11 X 17 

123C 

IBM 360/65 

Jordan 

0.00 

3.694 

0.532 


0.993 

2 X 15*^^ 

20 

UNIVAC 1106 

Widnall 

0.12 

4.073 

0.499 

0.444 


10 X 3 

73 

CDC 6700 

Cunningham 

0.12 

4.043 

0.524 

0.445 


6x5 

115 

IBM 370/155 

Cunningham 

0.12 

4.14S 

0.523 

0.443 


3x3 

7 

IBM 370/155 

Lopez 

0.00 

4.047 

0.550 


1.099 

16<*) 

20 

CDC 6600 

Notes: 









(a) (Not applicable) 








(b) Rubbert's 
wing at « 

program Is nonlinear. The above results were obtained by running the uncambered 
■ 0.1 radians. 

(c) Each of the Rubbert runs for the swept wing took 117 seconds. 



(d) Chordal modes, 
(t) Spanwlse modes. 








(f) Rounds off the kinjc at the centerline for this run. 




(g) CDC 6700 tllne Is 

In' term of time on a COC 6400. 





From Langan and Wang 

^Ref. 27 gives a value of 0.555 for PMCP. 
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■root chord, and bending moment center of pressure CBMCP) measured from the 
wing centerline divided by semispan are presented for each program. In Table V, 
the lift derivatives are given In place of the lift coefficients. 

Fig. 29 presents the spanwise lift distribution for the tapered wing 
at 11.4 degrees angle of attack. Fig. 30 shows the spanwise distribution of 
chordwise location of the center of pressure for the tapered wing at 11.4 
degrees angle of attack. 

Most of the results of the present method are In fair agreement with 
the computer results obtained by others and with the experimental results. The 
most noticeable disagreement is the spanwise lift distribution shown in 
Fig. 29 where the values of the present method are higher than the others towards 
the wing tip. A similar trend was noticed In the case of a rectangular wing 
reported In Case 3. 

Case 5; swept wing with control surface subson I c unsteady. - Case 5 
consisted of calculations on a swept wing with partial -span flaps (planform 
No. 5) shown in Fig. 10 which compare with the results reported by Stark (ref. 
24). The oscillating modes for wing and control surface and the corresponding 
aerodynamic coefficients are defined In Example 2 found In Appendix A. 

The breakdown of the wing and the control surface is shown in Fig. 31 
for the case of 40 col location points. Fig. 32 shows the arrangement of the 
swept wing with even strip width and Fig. 33 shows the arrangement with 
reduced sti ip width towards the tip. The aerodynamic coefficients calculated 
by using these arrangements and for different numbers of collocation points 
are compared with those of Stark (ref. 24) and with the results obtained by the 
doublet lattice methods (obtained by private communication from J. Griffin, 

Vought Aeronautics Co.) in Fig. 34 , and in Table VI. It is difficult 
to compare the relative computational efforts of the different methods used for 
the results reported here due to the lack of Information available. 

In general, the results for the swept wing are not In as good agreement 
with the referenced values as In the case of rectangular wing calculations. 

The lift coefficients are In good agreement except at high subsonic speed. 

At Mach number 0.9, the imaginary value of lift coefficient Is considerably 
smaller than the referenced value. The real parts of moment coefficients 

^ , 

/ 

/ 

/ 
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CONTROL SURFACE CHORD = ,712X LOCAL SEMI CHORD 
WING ASPECT RATIO = 2 



Fig, 31 Arrangement of Aerodynsunic Elements for a 

Swept Wing with Partial-Span Control Surface 
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Table VI. Lift and Moment Derivatives for Aspect Ratio 
Swept Wing with Partial-Span Flaps 
(Leading Edge Sweep Angle=60°, Taper Ratio=.2376) 
Rotating about Root Midchord (K = .5) 




Present Method ! 

Reference 

ffech 

Srm 

^RP 

Strip 






No. 

Width 

Real 

Imag. 

Ref. No. 

Real 

Imag. 


Qi2 

16 

Even 

1.9277 

2.7281 

Stark 24 

1.7399 

2.6312 

0 

II 

64 

II 

1.8742 

2.6369 

D.L. (A)* 

1.8039 

2.5424 


II 

64 

Reduced 

1.9179 

2.6833 

jfrrfe 

D.L.(B) 

1.8180 

2.5493 


Q22 

16 

Even 

.2700 

1.019 

Stark 24 

.2335 

1.1844 


II 

64 

II 

.2400 


D.L. (A)* 

.2987 

1.221 







** 




II 

64 

Reduced 

.234 

1.010 

D.L.(B) 

.3034 

1.2202 


^12 

14 

Even 

1.9394 

2.4229 

Stark 24 
* 

2.3309 

2.9015 


II 

16 

11 

2.6948 

2.8539 

D.L. (A) 

2.3620 

2.7519 


II 

40 

II 

2.5694 

2.2033 

frA 

D.L.(B) 

2.3869 

2.7426 


IT 

64 

M 

2.5962 

2.8341 





It 

64 

Reduced 

2.6566 

2.8936 





Qi3 

14 

Even 

.3846 


Stark 24 

.0984 

.7658 


II 

40 

11 

.1412 






'5i4 

14 

Even 


.3133 

Stark 24 

.8003 

.0693 


II 

40 

II 


.0211 




.7 

Q22 

14 

Even 

.2234 

1.3227 

Stark 24 

is 

.3495 

1.5228 


II 

16 

II 

.4401 

1.254 

D.L. (A) 

.4326 

1.5172 


II 





** 




40 

II 

.5655 

1.2922 

D.L.(B) 

.4484 

1.5313 


II 

64 

II 

.3857 

1.277 





11 

64 

Reduced 

.3738 

1.298 





^23 

14 

Even 

.1276 

.5684 

Stark 24 

.0518 

.5608 


II 

40 

II 

.0153 

.5641 





^24 

14 

Even 

.5332 

.3479 

Stark 24 

.5663 

.1917 


11 

40 

II 

.5693 

.1397 
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Table VI. (Cont’d) 



* D.L. (A): Doublet Lattice Calculations by Rodemich's Foimulation 
**D.L.(B): " " " " Landahl^s " 

(Results obtained by private connunication) 
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are In close agreement with the referenced values up to Mach number 0.9; 
imaginary parts are, however, consistently lower than the referenced values, 
giving the worst comparison at Mach number 0.9 as shown in Fig. 34. Improve- 
ments were noticed as a higher number of col location points was used at Mach 
number 0.9. 

The flap moment coefficients were calculated with only 14 collocation 
points for Mach numbers ranging 0 to .9. Considering the small number of 
collocation points used, the results are In fair agreement as Is shown In 
Fig. 35. One calculation was made with 40 collocation points for Mach number 
0.7 which resulted in marked Improvement over the results obtained from 14 
collocation points. This indicates convergence of the present results to the 
referenced values with a higher number of collocation points. 

Case 6; T-tails - subsonic unsteady .- Two T-tails found In the 
literature were re-analyzed, using the first order zero pressure rectangular 
elements. The first was the rectangular T^tail (planform No. 6) shown in 
Fig. II. Results are compared In Table VII with experimental results reported 
by Clevenson and Leadbetter (ref. 28) and with calculations by Kalman et al. 

(ref. 23). In the latter calculations, the effect of the tunnel, wall as a 
reflecting plane was disregarded. In comparison, Kalman^s calculations used 
100 collocation point which is equivalent to using 200 points if the effects of 
symmetries had not been included in the computer program. 

The second was Stark’s (ref. 29) we I I -documented T-tail (planform No. 7) 
shown in Fig. 12. Results are shown in Table VII! compared with results reported 
by Stark and results by Kalman (ref. 23), using the doublet lattice method. The 
numbers of col location points used for the present method were 18, 50, 60, 66, 

70, and 72. Stxfy-three collocation points were used In Stark’s configuration 
and 230 in Kalman's. These numbers would have been higher if effects of symmetries 
had not been considered.. 

The computed results general ly showed convergence towards the referenced 
values as the number of collocation points was increased up to 60 collocation 
points. Further increase of the number of col location points to 66 and 72 
initially gave erratic results. However, when the elements adjacent to the 
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Tabla VIII. Canpartson of Aerodynamic Coefficients for Stark's T-tall^* 
(Mach No. » .8) 





Yaw (j 

= 1 ) 

■■ ^ 

Sidesway (J=2) 

Rol 1 

(JO) 

1 Otrfvattves 

» 

K 

r 

Real 

Imag. 

Real 

Imag. 

Real 

Imag. 


Present, 18 

.2 

.021 1 

-.3919 

.0367 

-.0053 

.0098 

.021 1 


Present, Nj^=50 

,2 

-.2849 

-.4734 

.02713 

-.06099 

.01244 

.02705 


Present,Npjp=60 

.2 

-I .221 

-.4316 

.03197 

-.03147 

.01253 

.02613 

Yaw 1 ng 







Moment 

Present, Npp=66 

.2 

-.918 

-.5860 

.9731 

-.4333 

.1083 

.00026 

(C|j) 

Present, Npp=70 

.2 

-.2719 

-.4950 

.02956 

-.05899 

-.01339 

.02893 


Present, Npp=72 

.2 

-.7165 

2.1534 

-.4507 

.404 

-.0029 

.0350 


23 

Kalmon^^ 

.2 

-.0837 

-.5270 

.0470 

-.0278 

.0137 

.0257 


Stark^* 

.2 

-.0961 

-.481 1 

.412 

.0300 

.0125 

.0239 


Present, Npp= 18 

.2 

-.5861 

-.3520 

.0265 

-.1177 

.0133 

- .0298 


Present, Npp=50 

.2 

-.7818 

-.3496 

.01269 

-.150 

.01662 

-.0362 


Present, Npp=60 

.2 

-.6231 

-.3446 

.02160 

-.1228 

.01653 

-.034 7 

Side 








Force 

Present,Npp=66 

.2 

-.3766 

1 .0956 

.3863 

-.6203 

-.4012 

.1477 


Present, Npp=70 

.2 

-.7793 

-.3602 

.01468 

-.1495 

.01766 

-.03592 


Present, Npp=72 

.2 

-.7991 

1.2299 

-.2777 

-.0483 

.00833 

-.0297 


Ka 1 man^^ 

.2 

-.6270 

-.3965 

• .0297 

-.1260 

.0171 

-.0318 


Stark^' 

.2 

-.6108 

-.3625 

.0241 

-.121 1 

.0158 

-.0295 


Present,Npp=l8 

.2 

-.1262 

-.1101 

.0129 

-.0264 

.0120 

-.0382 


Present, Npp=50 

.2 

-.1231 

-.1346 

.0177 

-.0266 

.0201 

-.0596 

Rolling 

Present, Npp=60 

.2 

-.1295 

-.1354 

.01693 

-.02769 

.01967 

-.0582 

Moment 

Present, Npp=66 

.2 

-.4012 

.1477 

-.0465 

-.0632 

.0142 

-.0625 

(C3J, 

Present, Npp=70 

.2 

-.1236 

-.1373 

.01815 

-.02657 

.02065 

-.0594 


Present, Npp=72 

,2 

-.07963 

-.3133 

.0512 

-.0310 

.021 1 

-.0586 


, 23 

Kalman 

.2 

-.1270 

.1267 

.0154 

-.0269 

.0186 

-.0529 


Stark^^ 

,2 

-.1247 

-.1151 

.0134 

-.0255 

.0179 

-.0497 

Yawing 

Present, Npp= 1 8 

.3 

.0458 

-.6220 

.0844 

-.0160 

.0244 

.0331 

Moment 

29 

Stark^ 

.3 

-.0690 

-.7736 

.0973 

-.0549 

.0315 

.0403 

Side 

Present,Npp»l8 

.3 

-.6106 

-.5145 

.0552 

-.1818 

.0300 

-.0458 

Force 

29 

Stark^^ 

.3 

-.6471 

-.5592 

.0562 

-.1895 

.0379 

-.0449 

Rol 1 Ing 

Present.Npp«ia 

.3 

-.1342 

-.1613 

.0269 

-.0428 

.266 

-.0588 

Moment 

29 

Stark^"" 

.3 

-.1344 

-.1729 

.0299 

-.0415 

.041 1 

-.0770 


* BaMd on y b (Sm Fig. 12) 
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intersection of the two surfaces were kept to the same size as for the case 
of 60 collocation points, the results obtained with 70 collocation points 
showed a tendency to converge. 

The evaluation of the three versions of rectangular subsonic elements 
for planar cases showed that the best agreement resulted from the "first 
order, zero pressure, rectangular element" followed by the "first order, 
constant potential, rectangular elements," with the "zero order, constant 
potential, rectangular elements" appearing to have the least agreement. 
These three methods were applied to the T-tail with 18 elements, but the 
results given in Table IX do not indicate that any one method is preferable 
over the others. 

Definitions of aerodynamic coefficients used for the two T-tall 
calculations are as follows: 

F 

C , = 1 — 

y,ip irq 


2 

where = Side Force 


M 




Pitching Moment 


^FIN 


Fin Area 


Stark’s T-ta i I ; 


* j 


Tiq b 


;/ 

Tall 


J (x,y)AP (x,y)dxdy 


+ // 
Tall 


J (x,z)AP (x,z)dxdz 
m ' n ' 
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. 29 

Table IX. C3citparison of Various Aerodynamic Elements for Stark's 

T-tail (N_ = 18, K = .2*) 

RP ' r 


Mode 

Derivatiw-^^^,^^^ 

Yaw (j=l) 

Sidesway (j=2) 

Roll (j=3) 

Real 

Imag. 

Real 

liiHg. 

Real Imag. 

Yawing 

Moment 

(C„) 

By 

ZOCP 

0.321 

-.4222 

.0416 

-.0045 

.0114 

.0214 

By 

FOCP 

.0289 

-.4080 

.0389 

-.0044 

.0104 

.0215 

By 

FOZP 

.02094 

-.3978 

.0367 

-.0053 

.0098 

.0211 

By 

Stark 

-.0961 

0.4811 

.0412 

-.0300 

.0125 

.0239 

Side 

Force 

(C,.) 

By 

ZOCP 

-.5671 

-.3575 

.0289 

-.1141 

.0144 

-.0289 

By 

FOCP 

-.5787 

-.3548 

.0275 

-.1163 

.0138 

-.0294 

By 

FOZP 

-.5859 

-.3519 

.0265 

-.1176 

.0133 

-.0297 

By 

Stark 

-.6108 

-.3625 

.0241 

-.1211 

.0158 

-.0295 

Rolling 

Moment 

By 

ZOCP 

-.1240 

-.1083 

.0127 

-.0257 

.0122 

-.0377 

By 

FOCP 

-.1252 

-.1094 

.0128 

-.0261 

.0122 

-.0380 

By 

FOZP 

-.1260 

-.1100 

.0219 

0.0263 

.0120 

-.0382 

^ 31 

Stark 

-.1248 

-.1151 

.0134 

-.0255 

.0179 

-.0497 


* Based on ^ b (See Fig. 12) 

Effective wake length (F=6) 

ZOCP = Zero -order. Constant Potential, Rectangular, Element 
FOCP = First order. Constant Potential, Rectangular, Element 
FOZP = First order. Zero Pressure, Rectangular, Element 
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where the oscillating modes considered are 

n=I 

n=2 for tai I 

n=3 

, n=l 

, n=2 for fin 

, n=3 

I = I * ’ *m 
j = I **‘n 
c = b/3 

Xq = x-coordinate of fin tip mid chord 
AP^ = Pressure jump correspond i ng 

Case 7: rectangular wings - supersonic steady .- Using in-plane steady 

supersonic rectangular elements, two sets of calculations were made on rectan- 
gular wings of various aspect ratios. The first was for comparison with results 
previously obtained by Haviland (ref. 30) for Mach numbers ranging from 1/2 to 
[.8 and aspect ratios as shown in Fig. 13. The second set was made to 
compare with the results reported by Nelson idt al. (ref. 31) for Mach number 1.3 
and aspect ratios as indicated in Fig. 13. 

The breakdown of the wing surface for supersonic calculations is different 
from that for subsonic calculations. This has been illustrated in Fig. 36 for 
the case of 16 collocation points on an aspect ratio two rectangular wing. The 
lift and moment derivatives calculated for the rectangular wing pitching about 
the leading edge which compare with results by Haviland (ref. 30) are shown in 
Table X. The derivatives calculated to compare with the results reported by 
Nelson et a 1 . (ref, 31) for rectangular wings pitching about the midchord are 
shown i n Tab I e XI . 

Using 8x8 arrays of elements, the present results are within 5^ of. the 
referenced values for the whole Mach number range used in the calculations. 


J (x,y) = 
n ^ 


0 

0 

y/b 


J^(x,z) = 


(Xq - x)/c 


-z/b 
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Arrangement of Aerodynamic Elements 
for Aspect Ratio Two Rectangular Wing 
for steady Supersonic Calculations 


Fig 


36 


Method 


No. of 
Points 


Mach 

No. 


Havi I and 
Present 


30 


6x6 
4 x4 
5x5 
6x6 
8x8 


1.2 


Havi I and 
Present 


30 


6x6 

4x4 

5 x 5 

6x6 

8x8 


1.4 


Havi land 
Present 


30 


6 x6 
4x4 
5x5 
6x6 
8x8 


1.6 


d Moment Derivatives 
for Rectangular Wings 
Flows 


■ 


Aspect 

Ratio 

Lift 

Derivative 

Moment 

Derivative 

3.0151 

4.49 

00 

0 

• 

1 

II 

4.878 

-4.4056 

II 

4.8158 

-4.3398 

It 

4.7507 

-4.2595 

II 

4.694 

-4.1979 

2.0412 

3.08 

-2.77 

H 

3.2922 

-2.9683 

II 

3.2503 

-2.9246 

II 

3.2161 

-2.8836 

II 

3.1589 

-2.8159 

1.6013 

2.42 

-2.17 

11 

2.5827 

-2.3286 

II 

2.5499 

-2.2943 

It 

2.5230 

-2.2622 

II 

2.4984 

-2.2369 




















Table X 


Cont ' d 


Method 

No. 

Points 

Mach 

No. 

Aspect 

Ratio 

Lift 

Derivative 

Moment 

Derivative 

HavH and^^ 


1.8 

1. 3363 

2.01 

-1.79 

Present 


II 

II 

2.1555 

-1.538 

II 

5 x 5 

II 

II 

2.1279 

-1.9147 

It 

6 x 6 

It 

II 

2.1054 

-1. 8876 
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Table XI. Comparison of Lift and Moment Derivatives (about 
Midchord) for Rectangular Wings in Steady Super- 
sonic Flows Mach No. - 1.3 


Method 

No. of 
Po i nts 

Aspect 

Ratio 

Li ft 

Coeff icient 

Moment 
Coef f icient 

Nelson^* 
et a 1 . 


1 .2039 

2.42 

.40 

Present 

4x4 

tt 

2.596 

.465 

M 1 31 

Ne 1 son 





et al . 


2.0 

3.38 

.245 

Present 

4x4 

It 

3.678 

.234 

It 

8x8 

ft 

3.494 

.254 

Nelson^* 





et a 1 . 


4.0 

4.16 

.120 

Present 

4x4 

tt 

4.347 

.105 

IT 

8x8 

tt 

4.229 
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Definitions of the derivati vesL stiov^n in the tables are the same as those 
used in Case 2. 


Computation Time 

The four-by-four results took betv/een thirteen Cl 3) and twenty seven (27) 
seconds to run on the CDC 6400, using batch processing. The eight-by-eight 
results required 220 to 294 seconds. The time should be roughly proportional 
to the number of elements in the aerodynamic influence matrix, i.e.., to the 
number of reference points (N„„) squared, or 1:16 in the above case. 

Kr 

PROPOSED EXTENDED METHOD OF ANALYSIS 

In this section, the method described previously is extended to polygonal 
elements for subsonic and supersonic flows In the out-of-plane case. Expressions 
are shown to be the same as those previously developed when applied to rectangular 
e I ements . y 

Statement of Problem 

The problem of the extended method is identical to the original problem. 

Coordinate System 

This is Identical to the original system outlined. However, for analytical 
convenience, use Is made of the X,S,N coordinate system defined in Eg. (47). 

This gives the coordinate of the sending point in a system v/hose origin is the 
receiving point Jr. X is parallel to the x’ axis, S is parallel to the plane of 
sending surface, and N is normal to It. 

Governing Differential Equation 

This is identical to Eq. (2). 

Integral Equation 

This is defined by Eqs. C4) and (5). 
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Expression for Greenes function. - Ttits Is equivalent to Eq. C6). However, 
the cutoff function. is redefined in terms of a real argument, and the results 
are expressed in terms of the X,S,N coordinate system. 

Taking jr at the collocation point k, the subsonic Green’s Func'ion Is 


G = Gp{frCk),(P} = CR + MX]} 

Define also 

G^{lr(k)jp} = ^ exp C-R + MX]} 

where 

' R =a/x2 + 3^s2 + ^ 2^2 


The supersonic Green’s function may then be expressed as: 

G = U(F)CGp^{ir(k),fp} + G^ {ir(k),jp}] 

U is the step function, defined by 


U(x) = 


I; X > 0 
0; X < 0 


(79A) 


(79B) 


(80) 


and the cutoff function F (which is positive only in the forward Mach cone) is 

F = -X 

as used by Watkins and Berman (ref. 14). 

Expression for This is still defined by Eq. (5), but {t j§ expressed 
in a different form than in Eq. C7). Note first that 
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3/an = -n - v3/3v - n • cr8/3a 


^ ~ -cosCy - g)3/3v - sinCy - g)3/3(j (81) 

where y and are the angles between the local normal and the z-axis on the 
sending and ocefving surfaces respectively, as shov/n In Figure I. Then 

= ~cos(y “ g)3^G/3v^ - sln(y “ g)3^G/3v0cr (82) 

Method of Solution 

The method of solution is generally the same as suggested previously. 
However, expressions for irregular elements are developed for both subsonic 
and supersonic flow. No doubt the same can be done for sonic flows by deriving 
the limiting expressions correctly. 

-Integration by Discrete Elements 

The expression for the "aerodynamic influence matrix" ^^remains as 
in Eq. (12). 


Aerodynamic Elements 

The catalog of aerodynamic elements is extended. Due to the new form 
for expressed in Eq. (82), however, the new and old elements can only be 
compared for the in-plane case, 

Basic form .- The zero-pressure form of element has given the best results 
up to now, and has therefore been chosen as the fundamental form for the 
extended method. Also, Eqs. (13) and (34) have been combined, f has been 
redefined, and the X,S.,N, coordinate system has been introduced, resulting in 


'A<f(il*,X,S) = exp 



CX - X(Jl)D 


(83) 
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The net pressure differential AP acting in the positive N direction Is 
now given by 



j 3^ iw( A(j)0>*,X,S) 
( 3X V ) V 


(84) 


Therefore, for any element on which f equals I, AP Is zero, (the zero 
pressure element), while for any element on which f equals {x - x(A)} 


= -2^ expj^ CX - X(£,)D| (85) 

which is constant In the zero frequency case. The relative advantages of 
using polynomial forms for f are not understood at this time. However, the , 
above two examples warrant further discussion. 

In the first case with f equals 1 there is no pressure except at the edge 
of region A* where there is a lift line correspond ing to the stepwise jump In 
A(J). This would be treated as a vortex loop in the steady state case, as is 
done for example, in the method of Rubbert and Saaris (ref. 5). 

In the second case, the pressure varies harmonically. Further possibili- 
Itles are polynomials with unknown coefficients. For example. If f Is 
1 + a {x - x(5,)}, we have two unknowns for each element, A4>/v and a. 

The expressions for the general terms A. of the aerodynamic Influence 
matrix given in Eqs. (14) and (15) are retained. It might be noted that the 
form of Eq. (14) could be used for the on-surface elements as well as for 
the wake region elements if it were desired to subdivide them. Also, the form 
of Eq. (15) might be used for trailing edge elements, if the ^panelled' wake 
used in earlier parts of this report were to be replaced by integrals extending 
to infinity. 

It is thus convenient to introduce a modified Greenes function 


G = G exp {- iy [X - X(Z)J } 


( 86 ) 
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Then using Eqs. (82) and (83), Eq. (16) can be rewritten as 


k^(k, «.,£.*) - “//jj,* {cos(y - g)82G/aN2 

+ sln(Y - g) 32 G/ 3 N 3 S} f(«,*,X,S)dX dS (87) 

Subsonic aerodynamic elements .- After some manipulation. It can be shown 

that 


3^G _ r (, ^ |lj _ RG (88) 

~ I / ^ oTT < <!>mZ \ D { TTTdT I 

and 

fj, 

aN3S I I' 

where 

V G^exp I - - X(4)] I 


p-'C 


dIN' 




J 


Io)R I 3^ 

j 3N3S 




(89) 


Minimum first order expansion-subsonic ; The argument of the exponential 
In the expression for Gp, Eq. (90), can be expanded about the point 1* as 
before for the development of the so called "first order" elements. Since the 
argument of the expanded exponential Is of the form luA/c (where A Is a typical 
dimension of the element £*), expansion up to the first order term should 
be adequate when smal I elements are used. Then, from Eq, (90), 


Gr ^ ik i' - CR - R(k,**>: - ^ - x(^»)d 

I 


X exp ^ CMR(k,il*) + X(£*) - B^X(Jl)3 | 


(910- 
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By substituting this into Eqs. (87), (88), and (89), the '’Minimum First 
Order Expansion" results in which ai i terms in and higher order are dropped* 
This leads to the following expression for k 


k. (k 


I Si*) = [ | cos(Y-g) j + s l n(Y-g) I t»(k £*) I 

' I ) 4ir o 4ir o j 

I + + ia)X(«,*) I _ ( cos(Y-g) 

^ j j 4 tt ^ 


sin(Y-g) T»(k 9,*)! |4^! 

47T 


exp |- CMR(k,£*) 


+ XU*) - e^xCA) 


]| 


(92) 


In the above expansion, terms such as (ojN/c)^ are dropped. The 
validity of this step is in question and is discussed later. The remaining 
terms dropped are of the form (ojA/V )2 which are small if the elements are 
themse I ves sma I I . 

The "I" integrals have the form 


= // f(**,X,S) ^ jij 

z* 


dXdS 


(93) 


Ij'tk.Jl*) = II f«*,X,S>X ^ 4j 
i* I ) 


dXdS 


(94) 


I^"(k,il*) = // f(Jl,*X,S) 

Z* 


9N9S (Rj 


dXdS 


(95) 


= II f(t»,X,S)X 3^ jij dXdS 
z* ^ ^ 


(96) 
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S upersonic aerodynamic elements ^- As wl+h Eqs. (88) and (89), It can be 
shown that 


^ ""^18 


8N3S 


•» 


1 


i 


■ 

1 I io)R j 

< 

(' -6^1 


32 

w 


loj S^u 


o)^UN2 




B^c sN- 

io) 3^0 oi^UNS 

B2c 3N3S " c2r3 


RG|R 


RG R 
A 


( 97 ) 


(98) 


iR 


The notation- [/^ refers to alternate choices of the subscripts, R (retarded) 
and A (advanced). The corresponding signs on the RH side are indicated by 
± or + as appropriate. 

M inimum first order expansion-supersonic : Using the same procedure as 

in the subsonic case, an expression is obtained for which corresponds to 
that given by Eq. (92). 


k.(k,£,JL*) = 


I 4 tt o' 4 tt o 


.. coR(k,A*) , 2ajR(k,ii*) o)R(k,£*) ^ 2i(oX(£) o)R(k,Jl*) 

X (2cOS -9 + 7T sin -A + r cos ^ 

I b2c B^c b2v g2c 


i cos(Y-g) T f 

I 47T 


+ sin^pli;'(k,J,^) l^cos 


2‘uii ^ o)R(kj^.*) 




j 

1 


J J"(k,^*)} sin — 


4ir 




B-^c 


lu 


X exp <- LXU*) - B^XU)2 J + k^”(k,£*) 


(99) 


The first term includes all contributions except from coincidence of 

the point k with the element £.* which is covered by k.'*. The derivation of 

<P 

the corresponding integrals has been discussed. The second term which covers 
the coincidence effect appears only in the supersonic case. 
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Although this equation appears formidable at first sight, it has a 
symmetry which should make programming relatively simple. 

The integrals are 


<(k.z*> = // fa*.x,s) ^ j-yp-j 


dXdS 


Il'Ck.t*) = // jy^l dXdS 


z* 


= // f(4*,x,s) ^ te 

z* ' 


dXdS 


Ii"(k,£*) = // f(£*,X,S) X 


a2 


dNdS 


|T| 


dXdS 


J'(k,£*) = // f(£*,X,S) dX,d$ 

3N2 


J’’(k,£*) = // fC£*,X,S) dXdS 


£* 


(100) 


(101) 


(102) 


( 103 ) 


( 104 ) 


( 105 ) 


The first four ”1" integrals differ from the corresponding subsonic ones 
only in that they contain the step function U(F), The "J" integrals only appear 
in the supersonic case and their integrands are zero except on the forward 
Mach cone. 


Evaluation of Integrals by Contour Integration 

Evaluation of the "I" and ”J" integrals begins with their reduction to 
contour form. These contours will be assumed to be polygons, so that it will 
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be necessary only to evaluate the contributions of straight line segments. 
Define the indefinite F integrals as fol lows: 

= // XV ^ j dXdS (106) 

<I07) 

It will be shown that these can be reduced to line integrals, thus, for 
example, if f (A*,X, S> equa I s unity, then Ij can be expressed as the sum of 
line integrals along the sides of the polygons, so that 

N ( „ ) I+l 

Il"(k,«,*) = J , ( 108 ) 


where i, i + I are two consecutive points (in the right handed sense) around 

the contour, as shown in Figure 37. A transformation to an X,’ S,' N 

coordinate system will be used to facilitate integration. The bracket 

{F}i"*^* has a particular sighif icance, If the line i to i + I does hot cross 
^ i + 1 

the Mach cone. It is identical to LFJj , and has the conventional meaning 
"evaluate F at i + I and subtact value at i." However, there is no contri- 
bution to the line integral from outside the Mach cone, but there is one from 
the intersection point, so that a different Interpretation must be made when a 
line intersects a Mach cone. 

In Eq. (108) above, the shape factor f corresponds to the zero pressure 
case which has already been discussed. In the harmonic pressure case, with 
f equal to x - xOl), one would have 





X(j,*)F 


If 

1,0 


i+l 

i 


(109) 


Given a method to evaluate the "P” integrals., a large range of elements is 
available. Five cases of Interest can be identified concerning two consec- 
utive points, I and 2 which are as follows: 
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Contour Integral 


(t) Subsonic case, and supersonic case with points I, 2 lying inside the 
Mach forward cone. 

(II) Supersonic case with point I Inside the forward Mach cone and point 
2 outside. 

(Ill) Supersonic case with points I and 2 reversed relative to I I . 

(IV) Supersonic case with points I and 2 outside the forward Mach cone, 

(V) Supersonic case with collocation point k lying on element I*. 

Cases I and II are discussed in the following subsections. Case ill 
is easily deduced from II. Case IV yields zeros (for this line section, though 
not necessarily for the whole element). Finally, Case V Introduces additional 
diagonal terms. (Note that a point lies inside the forward Mach cone If the 
cutoff function F defined under Eq. (80) is positive). 

The following transformation of coordinates is used in the plane 
element I* 


X = “ Xi) - 32 SMS 2 - Sj)} /Ri^2 (1 10) 

S = Ej^2 {x’(S 2 - Si) + S'(X2 - Xi)} /Ri^2 (N!) 

dXdS = dX'dS' ( I 12) 


where Xi, Si, and X 2 , S 2 are the coordinates of the points I, 2 and 3 


Ri, 2 = V1(X2 - Xi)2 + ^2(32 - Si)^l (1I3A) 

and 


^1,2 " ' ^^2 ’ ^^^^2 ” ^ ° 

- -I If (X 2 - Xj)2 + 62(S2 - Sj)2 < 0 (II3B) 
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In the supersonic case, Ej 2 E!q, (113) is negative unity if the line 
from point no, I to point no, 2 makes a smaller angle with the S-axis than does 
the Mach line. Also, Ri ^2 is zero when the line between the two points is 
parallel to a Mach line. This case is discussed later. 

The Inverse transformat ion to Eqs. (NO) to (III) is 

X' = -|X(X2 - Xi) + e^S(S2 - Sj)} /Ri,2 (114) 

S’ = {-X(S2 - Si) + S(X2 - Xi)} /Ri,2 (115) 

In terms of the new coordinates, 

R= VEi,2X’^ + 62Ei,2S’2 + (116) 

In the subsonic case the transformation can be represented by the following 
three transformations app.l led in order: 

(1) A Prandt I -Gauert Transformation 

(2) A rotation of axes so that the line from points I to 2 is parallel 
to the X' axis. 

(3) An inverse Prandt I -Gl auert transformation. 

In the supersonic case, the correspond ing transformation which transforms 
the equations to rectangular hyperbolic is used (See, for example, the report 
by Heaslet and Lomax (ref. 32)), but the factor E^ ^ is Introduced to avoid the 
use of imaginary numbers. * 

It is read.! I y shown that 

dX = . - /.Xllox' - -- S -ii (117) 

f*l,2 R)^2 
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( 118 ) 


dS = ^l,2<g 2 - S)) 
Ri»2 


dX' + ^1 . . 2 <X 2 - X j. ) , jj, 
Ri,2 


also that 


I dX' = Ei,2Ri, 2 
J dS' = 0 


(119) 

( 120) 


Note that It Is possible for X* to be positive in the forward Mach cone, 
whereas X is always negative, as Illustrated in Figure 38. A useful method of 
evaluating Integrals can be derived from the following: 


{ 


3 ^ 




/47t 6(0, 0,0); M < I 
\ 2tt 6(0, 0,0); M > I 


( 121 ) 


given 


Case I, no Mach cone intersection ." 

”F'* integrals of the first kind ; These Integrals have the general form 
In Eq. (106), but with U(F) = i. 

To cover the zero pressure case, with f = I, we require m = 0, I; n = 0 


'^ 0.0 = //|^ 

Substituting from Eq. (121), using Green’s Theorem . (ref .21 ), and then Eqs. (114) 
and Cl 15) 


-// 1 

~ j .dXdS 


}w{ ^ 

fr jt] 

dX 

j^] " 


\k\ 


(122) 
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Mach cone 


Figure 38. Case of Positive X* in Forward Mach Cone 
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Taking the straight line from point no. I to point no. 2 which is 
parallel to the X* axis, dS* is zero. Thus, 



li 

3S* |R 



R(E 


i + N2)J 


2 

1 


(123) 


It is convenient to leave the Integral In this form since, for computer 
evaluation, X,* S* can be found readily. 

K.0 = // >^|^j 


Applying the same substitutions as before, together with Eq. (NO), and using 
the rule for differentiation of products, one has 


If' = 




Rl 2 


. g^Ei.^SMS^ - Si) ^ j^j ^ S^ .,2 (S ; - S j ) J 

^ f • 


3S’ R 


Ri,2 


dX* 
1 R 


3^S 

RR 


jEi.2(X2 - X.) 

^ e l »fr:r^^r .- S x> log X- + r[] ^ 


( 124) 


Note that the last term In Eq. (124) is complex if Ei,2 negative. 
Evaluation of the integrals in Eqs. (123) and (124) is facilitated by the 
standard Integrals of Appendix B. 


”F” Integrals of the second kind : These "F” integrals have the general 
form given in Eq. (107): 


" a2 

''o.o = // 


3N3S 


1^1 


dXdS 
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Integrating first with respect to dS 


^ W l^l '‘X 


then transforming and taking the straight line element from point I to 2 

2 ^ 2 


I . Ei,2(Xa - X;) _i 111 <jx' 

- ■ Ri,2 ^ 3N |R( 


Ei, 2NX»(^- Xi) 


2 


RR 




(b) 


1.0 


= //X 


82 (| 

8N8S |r 


J 1 

I 

dXdS 


As before, for the straight line element from point I to 2: 

2 


L" I El 2<X2 - f2 X' ^ lijdX’ 

(Rj‘2> j SNjRp^ 

, b^sMXz - Xi)(S 2 - Sj) a li ) 

•'l 3N 


(125) 


(126) 


^ r 62 n(x,> Xi) L . X^S*(So - $i) l] ^ 

L j^l,2^^2 (Ej^aF^ + |![2) JJ ^ 


(127) 


Case II with Mach cone intersection .- Let point C lie on the intersection 
of the forward Mach cone with a straight I ine path from points I to 2 so that 
point no. 2 lies outside of the cone. The integrands of the ”1" and '*J'* integrals 
are now zero beyond point C, but most of their terms are singular at C. This 
arrangement Is illustrated in Figure 39. 

First note, that when point no. 2 is outside the Mach cone, function U(F) 
goes from unity to zero at a point C which I ies between points I and 2 as shown 
-in Ffgure 39. Then, since the derivative of the step function Is the Dirac 
delta function 6(X), it follows that for any function Y(X) 
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Figure 39; Intersection with Forward Mach Cone 
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Jy(X) 3U(X)/aX dX = /Y(X) 6(X) dX = Y(0) 


In the more general form, applicable to the case in point 
/ Y(X,S,N) ~U {fCX,S,N)( dX 

= 'c (128) 

where a is X, S or N. Equation (116) for R can be factored in the form 


R =V-Ei,2F'F'a C29) 

where 

F'.= -X’ -V-3^S'2 _ ^n2‘ (I30A) 

F^’ = X’ - V-3^S’2 - (I30B) 


It is not easy to identify the cut-off function for the forward Mach cone. 
For the following derivation, it will be assumed that It is F’ as given by 
Eq. (I30A), which would be the case if points I and 2 were aligned in a positive 
direction along the X axis. 

Note that, at a point C on the forward Mach cone 


x' = -1 /-b2S'2 - s 2 n 2 
eye e 


First, the following indefinite integral is evaluated at C 


(131) 


/ 


3 

3S' 



/ 


|l 3U(F’) 
(R 3S’ 


32u(F’ )S' 



I 3^XVS^__ 

(r\ 43^S’2 - '32(^2 R(g2s’2 + ^ 2 .^ 2 .) 
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(132) 


e^Rs» 

X'(3^S’2 + 32n2 


= 0 


By in+erchanging variables 


3N 




= 0 


(133) 


Also, the fo! lovnng indefinite integral 


f X' -j- dX' I 

j ^ 3S' I R j l( 


3^5 ’X* 




3^U(F»)X'S' 


R/32S.2 _ 32^2 R3 


dXi 


lc = o 


(134) 


Again, by interchanging variables 


/X’ W 1^1 «X’ 1^ = 0 '135, 

It follows from Eqs. (132) to (135) that all of the singular terms in the "F” 
Integrals vanish on the Mach cone. Presumably, the same results would follow 
from the application of the method of determining the principal values of 
Improper integrals given by Mangier (ref. 33). 

*'F” integrals of the first kind with Mach cone intersection; 


(a) 


f 

Ro,0 


cr 9^ (UCFM 
'* ) R 


dXdS 


Proceeding as for Eq. (123) but substituting from Eq. (132) 




S ucfm) 


dX’ 


■ El ,.z X . ^S L. _ 
R(Ej^2S’^ + N^) 


(136) 
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Thus an element in the form of a polygon would have contributions from the 
'VeVtlces inside the forward Mach cone only. 


(b) 


cr' _ rr X, 3^ 1 U(F) 

^1.0 ^ ) R 


dXdS 


As for Eq. (124) but substituting from Eqs. (133) and (134) 

t )2 




1^. V/ ^ , X'S'CSs - Si) ) 

' '^’^1,2 {Ei^ 2S'“ + N^)j 

i^g {/eTJx' + r} j 


(137) 


Thus there is a contribution at point C from the finite term. 

integrals of the second kind with Mach cone intersection; 


(a) 


0 0 


= // 


8^ 

8N8S 


U(F) ) 


R 


dXdS 


AS for Eq. (126), but substituting from Eq. (133) 

.2 


P It ) _ El .?NX'(X7 - Xt) 

0,0 ' RRi,2^Ei,2S'2 + N2) 


(138) 


(b) 


1.0 


= //X 


U(F) 


8N8S 


dXdS 


Again, as for Eq. 


(127), but substituting from Eqs. 


(133) and (134) 



3^N(X9 - Xi) 
R(Rl ,2)2 



X^) + 


X’S'iSy - Si) 1 . 

(Ei^aS'^ + N^yj 1 


(139) 


”K" Inte grals of first and second kind: In order to evaluate the "J" 

IrtTegrals, v/e first define the "K" integrals which are similar to the "F” 
Integrals of Eqs. (106) and (107) 


■ loo ~ 


, ^ II xi"sn 32 u(F) 


m,n 


9N^ 


dXdS 


K ^ ff X^S 
m,n 


n 82U(F) 

8Nas 


dXdS 


Then, for example, if f(A*,X,S) = I , we have 





( 140 ) ' 

( 141 ) 


( 142 ) 


There is evidently no contribution except at a Mach cone intersection. 


(a) = // dXdS 


As for ^ in Eq. (123) 
0.0 




-3^S' 




3 2$ 

IT 


(143) 


Note that the sign would be reversed if the point I were outside the Mach cone, 
and point 2 were inside. 


(b) 


rr 32U(F) 
Kn.D = // 


3N3S 


As for Fq^o in Eq. (126) 


K 


0.0 


= Ei. 2(X, - Xi) ItKFll 
R. 3N 


1^2 


g2(X9 - Xi)N 


( 144 ) 
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Additional Supersonic Terms 


In evaluating the "I" and "J" integrals of Eqs. (100) to (105), the RH 
term of Eq. (121) was omitted. The term k^' in Eq. (99) replaces this, and must 
be evaluated for supersonic flows. 

In the original calculations for rectangular elements, the collocation 
points were placed directly over the lifting lines, as shown in Fig. 5. 
However, the procedure for calculating k" could be described as follows: 

The wing is divided up into rectangular winglets, as in Hedman’s method 
for the subsonic case, (see Fig. 14) but, in place of the lift lines at 
the quarter chords and the collocation points at the three-quarter chords 
of the winglets, the lift line and collocation point are both placed at 
mid-chord. Then, for the purpose of calculating the effect of a lift-Ifne 
directly on the collocation point with which it coincides, the lift line 
is replaced by the same total lift uniformly distributed over the element. 
Because we are actually concerned with an aerodynamic matrix v/hich relates 
velocity poten'. ial to downwash, the resulting matrix has some off-diagonal 
terms. The lift and downwash on a two-dimensional supersonic wing are 
directly related through the so-called "Ackeret result" (reference 4), 
so that the wing section consistent with a lift- line is stepped, nevertheless, 
the procedure outlined above gives correct lift and moment when applied to 
two-dimensional wings. As with the subsonic case, the simplified method 
used for rectangular elements will no longer be as effective when applied 
to arbitrary elements. 


Numerical Problems In Computation 

In evaluating the expressions for the "F" and "K" integrals, numerical 
problems may arise in the computer due to division by zero in such cases as 
the f o 1 1 ow i ng : 


(a) R = 0 

(b) S’ = N = 0 

(c) = 0 

It will be necessary to take care of these cases either by taking special 
precautions during preparation of input by programming additional computer logic. 

(a) R = 0. This occurs when there is a collocation point k on a vertex of an 
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element. This situation should be avoided. In the supersonic case. It may be 
that R Is very nearly zero because a vertex Is nearly on the Mach cone. If 
R < e, where e Is an arbitrary small number, calculations should proceed as 
though the Mach cone Is Intercepted. 

(b)S*=N=0. This occurs when the collocation point is on the extension of 
a side of an element. It can easily be shown from Eq. (115) that. If S* is 
zero anywhere on a straight line element, it Is zero all along It. If N - 0 
also, this could lead to division by zero in ail four integrals. In all such 

cases, however, the limiting contributions are zero. Therefore, on the computer 
terms with in the denominator can be disregarded when S’^ + small 


(c) 2 “ occurs when the line segment is parallel to a Mach line. 

Since Ri^2 appears in the denominators of X’ and S,’ the limiting expressions 
must be evaluated carefully. Resulting limits on Eqs. (123), (124), (126) and 
(127) are. 


where X' and S’ 
unity. 


’ 2 

*^ 0 , 0^1 


t^I o>| - 


Ls'rJ I 

W: 

[ NX’(X9 " Xi) 1 
R(S’)2 J 

[ - e^NS(x? - Xi 
RS’ 




(145) 


(146) 


( 147) 


( 148) 


are evaluated as for X’ and S’ but with Ri,2 taken as equal to 


Comparison of Extended Method with Previous Results 

As was pointed out earlier. It is only possible to make a comparison of 
the final expressions for the in-plane case. However, for parallel airifoTts 
(where angles 9 are equal), Eqs. (82) and (88) can be written as follows: 
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K , = •{ g^/R^ + lo)/cR f G 


II 


+ {-NVr 2 > i3gVR2 + 3la)g2/cR - o)2/c2 } G 


(149) 


This agrees with Eq. (7), after substitution from Eqs. (17) and (18) and taking 
equal to unity and equal to -N^/R. 

In the case of an in-plane airfoil, Eq. (92) can be written in a form 
directly comparable with Eq. (35) 


+ XU*) - g2x()l)j I 

Now from Eqs. (26) and (93) with f equal to unity 
Ij^(k,H*) _ l_ fr 3^ (\ \ 

— 7' jj 


( 150) 



//{^}dXdS 


= kgg(k,Jl*) 


(151) 


also, from Eq. (32) and (94) with f equal to unity 


^ lit ^\r/ 


4¥g2 


dXdS 
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( 152) 


Si* 


Substituting back from Eqs. (151) and (152) Into (150), it is found that the 
expression is identical to that given by Eqs. (19) and (35). Thus the 
extended method leads to expressions which are identical to those developed 
earlier when parallel, zero pressure elements are used. As a further verifi- 
cation, Eqs. (151) and (152) can be compared directly for the rectangular 
element shown in Figure 40. First, we evaluate according to Eqs, (108) and 
(123). The contribution of each point comes from two terms. Transforming, 
using Eqs. (114) and (M5), and noting that (X 0 - X/\) and (Sq - Sq) are both 
zero, we have the following: 


l^ik,Z*) ^ i_ 

4tt 4irR/\ 





4ttRaX^Sa 


4ttXaSa 


47tXbSb 


Rc 

47tXcSc 


R. 


D 


47TXnS 


D^D 


(153) 


where Ra is R evaluated at the point A. The above expression is identical to 
that for kgg given in Eq. (26). Using the same approach in dealing with 
and substituting from Eq. (124), the contributions of the first two terms cancel 
at each point. This leaves only the contributions of the logarithmic expression 
along the lines A to B and C to D. Thus, 

log(eS + R)1 - |og(-gS . R)j (154) 

-"a ^ 

It can, however, be shown that 


-log(“3S + R) - log( 3 S + R) 


therefore the above expression is identical to that given for k^j^ in Eq. (32). 
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Only the steady state In-plane case has been previously developed for 
supersonic flow. From Eqs. (46) and (99), we now have 

k^(k,£<^) = 2ks3(k,Jl*) = 2lQ»(k,A*)/4TT (155) 

Comparing Eqs. (153) and (56) It Is seen that the extended method Is compatible 
In the supersonic case also. 

In summary. It has been shown that the methods given In this section 
are Identical to those given previously when applied to In-plane, first order, 
zero pressure, rectangular elements In the subsonic case and to in-plane, 
steady state, rectangular elements In the supersonic case, (imputations based 
on both of these cases have already been evaluated. The extended method would 
give Identical results in the cases discussed. 

Computer Program 

The expressions developed here could be Included as new elements In the 
computer program described in Appendix A. However, It would be considerably 
less time consuming to precompute as much as possible and then to calculate the 
aerodynamic matrix from the precomputed data for each value of cd. An examina- 
tion of Eqs. (92) and (99) shows that the aerodynamic matrix Aj^ ^ could be 
written as fol lows: 

^ Io)L(k,A) 

+ k’^ (k,il,Jl) + J|(A*)k^(k,Jl,£*) (156) 

In the above expressions, five real matrices of the same size as A^^ ^ are 
required. The Indicated summation Is over wake elements which can be rectan- 
gular, and can therefore be evaluated more rapidly than polygonal elements. 

In the present program, the normal flow (down wash) at each collocation 
point Is assumed to be known, and generalized forces are calculated after the 
velocity potentials have been solved. The same could be done with the extended 
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method; however, since the aerodynamic elements can now be made to conform with 
structural elements, it might be better to determine location and downwash 
of the collocation point from deflections at structural nodes and, finally 
to compute equivalent forces at the nodal points. If zero pressure elements 
are used, the forces act along the perimeter. The force per unit length on a 
line element when taken right-handed is proportional to its projection along 
the a axis. It Is thus in the negative direction (i.e. opposite to the positive 
normal coordinate v) when the projection of the line is positive. 

Thus if dF is an element of force on the perimeter, as shown in Figure 41, 


dF = pv^Atj^exp 



EC - I da 


(157) 


Missing Term 

A first order expression is used to replace the exponential terms before 
integration. However, a term in (coN/c),^ where N is the normal distance from 
collocation point to sending surface, is ignored (this cannot be expressed by a 
contour integral). If it should later prove to be important, it may be neces- 
sary to evaluate the surface integral by its midpoint value using the approach 
described previously. 


Summary of Extended Method 

The principal aim of this section has been to develop the expressions 
for arbitrary out-of-piane polygon aerodynamic elements. Each expression is 
in the form of a contour integral so that it is only necessary to compute the 
contribution of each straight-line element. The expressions developed are 
generally applicable for subsonic and supersonic flow, but special steps are 
necessary when a line element intercepts the forward Mach cone from a receiving 
point. 

This method has been proven to be identical to the one developed and 
checked out earlier in this report when applied to in-plane surfaces. Thus, 
much of the experience acquired to date remains valid. 


108 





The problem of defining the Ackeret lift term is not discussed in detail 
because it is believed to require considerable numerical experimentation and 
careful programming. Certainly, the idea that the pressure at a given point 
contributes directly to the flow at that point is straightforward. Also, the 
problem of finding optimum ways of locating collocation points has not been 
discussed. In the applications to rectangular elements discussed earlier, 
special schemes were used for locating collocation points and the rectangles 
themselves. However, if the aerodynamic elements are to be compatible with 
structural elements, such favorable arrangements will no longer be possible. 
Computed results may thus turn out to be unreliable because the normal flow 
distribution corresponding to constant potential or zero pressure elements is 
singular on the edges of the elements on which there are finite vortices. If 
selection of collocation points were arbitrary, then it would be necessary to 
use more uniform potential distributions so that finite vortices could be 
eliminated. The polynomial forms envisioned in Eqs. (106) and (107) might 
possibly prove adequate for this purpose. 

CONCLUDING REMARKS 

In the present investigation, a numerical method has been developed on the 
basis of the linearized theory by which the generalized airforce matrix can 
be calculated for steady or oscillating wings In subsonic and supersonic flows. 
The use of the downwash-velocity potential relationship in preference to more 
commonly employed downwash-pressure relationship gives a somewhat simplified 
form of the kernel function. The approximations made In the evaluation of the 
integrals Involved avoid some of the difficulty incurred in evaluating the 
improper integrals. Thus, one of the advantages of the method proposed here 
is numerical simplicity in comparison to other existing methods. The present 
method has been tested quite extensively for rectangular elements in steady 
and unsteady cases in subsonic flows. Its application in supersonic flow has 
been limited to steady planar rectangular wings. The results of these calcula- 
tions are summarized in the following paragraphs. 

It has been demonstrated that accurate results can be obtained for cases 
of steady and oscillating planar rectangular wings in subsonic flows using 
the downwash-velocity potential method. Among the three types of aerodynamic 



elements evaluated, the beat results were obtained with the so called "first 
order, zero pressure" rectangular elements, followed by "first order, constant 

i 

potential" rectangular elements. The "zero order, constant potential rectan- 
gular elements, although considerably simpler than the other two types of 
elements provided the least favorable results. Results appear to be convergent 
with respect to the total number of boxes used. The wake can safely be termi- 
nated In five chord lengths or leas. 

The lift distributions on rectangular and tapered wings in subsonic flows 
Indicate that the use of rectangular elements gives excess lift towards the 
tips. It has been found that the distribution improves with shorter spanwise 
divisions towards the tips (given a limitation on the total number of elements 
to be used). All other results, including the distribution of the center of 
pressure, and various coefficients, have shown that this method gives results 
comparable to those resulting from other methods of analysis. 

Calculations for oscillating swept wings and for swept wings with flaps 
using rectangular elements have been compared with those made by other methods. 
It appears that, when the rectangular elements are used, the agreement of the 
swept wing results with the referenced values is not as good as in the case 
of rectangular wings. Improvement was noticed as smaller elements were used 
towards the tip, but satisfactory agreement has still not been obtained at 
high subsonic speeds. It is believed that these results indicate a need for 
elements conforming to the planform such as those that have been developed in 
the so-called extended method. The flap moment coefficients show fair agreement 
with the referenced values. However, the results are limited to 18 collocation 
points, except in one case where 40 points were used. 

Fair agreement has been obtained with other results published in the 
literature for two-T-tail configurations using midpoint integration for out- 
of-plane terms. As the number of elements is increased, however, agreement 
first improves and then gets worse. This appears to result from the use of 
mid-point integration for out-of -plane elements and indicates a need for the 
extended method. 



The results obtained for steady rectangular wings using In-plane steady 
supersonic rectangular elements are good, but convergence with the referenced 
values seems to be slower than in the corresponding subsonic cases. 

Because all of these calculations were made using a computer program which 
contained as subroutines all of the rectangular aerodynamic elements cited in 
this report, it is believed that the effectiveness of such a concept of 
aerodynamic elements was demonstrated. In addition, the objective of developing 
a method valid for subsonic and supersonic flows was partially accomplished. 

As is shown in the section covering the extended method, this objective has 
been rea I i zed conceptua I I y . 

The proposed extended method is capable of reproducing all of the in-plane 
calculations performed to date. It also appears to have the features necessary 
for better results on swept wings and out-of-plane cases (such as T-tails). 

It Is therefore concluded that the extended method Is potentially capable of 
handling all calculations relating to oscillating wings. 

It is believed that the proposed extended method will contribute greatly 
to the problem of handling flutter and to other related problems in automated 
structural design. In particular, the following three points are of signi- 
ficance: 

1. The computation procedure is based on the concept of aerodynamic 
elements. This procedure treats arbitrary panels on the aero- 
dynamic surface which can, therefore, conform to the structural 
surface panels. Thus, not only can the computed forces employed act 
at the structural nodes, but, also, the preparation of aerodynamic 
input data is simplified since it is so similar to the structural 
Input data. 

2. This method applies to a broad Mach number range so that one computer 
program can handle subsonic and supersonic flows. 

3. The calculations can be arranged so that If five matrices are 
computed first, then the flutter forces for the complete reduced 
frequency range can be calculated rapidly. 
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RECOMMENDATIONS 


It is recommended that the proposed extended method be programmed and 
hested. Fol lov/ing this, the tasks listed below should be performed, 

(1) Adaptation of the method to the unsteady transonic (M = i ) caset 
Although this adaptation may not result in accurate air force 
predictions and its usefulness may thus be limited, it is needed 
to provide a basic linear solution capability against which more 
refined methods can be compared. 

(2) Investigation to determine the best method of including the 
Ackeret or "two-dimensional" lift effects in the supersonic case. 
The method used for the numerical results given in Tables 10 and tl 
was selected for convenience but is too restrictive for use in the 
aerodynamic element concept. 

(3) Determination of optimum chordwise spacing of elements. Thus far, 
the Hedman (ref. 2) method of locating lift lines has been used 
exclusively In calculations. Again, this is too restrictive for 
application to the aerodynamic element concept. 

(4) Development of a set of rules for selecting aerodynamic element 
configurations to meet given accuracy objectives. To date, only 
a few variations have been made on rectangular wings and those 
were made mainly to observe gross effects and to demonstrate 
convergence. 

(5) Trial inclusion of the method in an automated design computer 
system. One of the goals in developing this method was to make It 
possible to compute aerodynamic forces at structural nodes. There- 
fore, an evaluation of its capability in this application and of 
its practical application in general would be of great interest. 

(6) Trial development of new aerodynamic elements. It is especially 
important to investigate arrangements which will eliminate finite 
vortices since it is expected that they will remove restrictions 
on collocation point locations. 
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The Computer Program 

Introduction .- In order to perform the calculations on the downwash- 
velocity potential method, a computer program has been written In FORTRAN IV 
language to be used with CDC 6400 at the University of Virginia Computer 
Center, its principal feature Is a subroutine to calculate and store the 
aerodynamic Influence matrix which Is then Inverted and used to calculate 
generalized forces which result from specified motions of the surfaces. 

The program assumes that the user will supply data In the form of a 
breakdown of each surface Into a large number of aerodynamic elements. To 
date, eight sets of subsonic compressible aerodynamic elements (all rectan- 
gular) have been formulated. For uniform size rectangular elements, the 
required Input data can be computed by using subprograms from basic Informa- 
tion such as the geometry of the surfaces, their planform arrangements, and 
the specified motions of these surfaces. 

Program Descrl ptlon .- 

Summary of Programs; 

Main Program : 

FLUTT 

Subprograms : 

DMATR and DREAD 
WMAT and WREAD 

ACALC and A(X)N3 through AC0N8 and AGON 1 1,12 

CMATR and CC0N5, 6 and CC0NII,12 

ALIST 

TEST 

INVER 

Their roles are the following: 

FLUTT: 

•Reads In data controlling the desired output, controls selection 
of the subroutines to be used, and specifies Integer words, floating 
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words, and the I Ine counts for each type of data. 

•Calculates the addresses for various matrices. 

•Reads in Mach numbers and frequencies. 

DMATR: 

•Generates D vector for uniform size rectangular elements In each 
strip from the basic information. 

DREAD: 

• Controls reading in D vector for rectangular elements of nonun I form 
size. 

WMAT: 

•Generates W matrix from the modal displacements of each strip. 
WREAD: 

•Controls reading in W matrix for rectangular elements of nonuniform 
size. 

ACALC; 

•Controls retrieval of Data vector using subprogram TEST. 

•Determines the conditions of symmetry for each surface. 

•Calls for appropriate ACON’s to be used. 

•Generates A (aerodynamic) Matrix. 

AC0N5 through AC0N8 and AC0NII,t2 

Calculates terms In A matrix by various methods described 
previously. 

CMATR 

•Controls retrieval of Data vector by using subprogram TEST. 

•Controls retrieval of A and W matrices. 

•Calls for appropriate CCON’s to be used. 

•Calculates velocity potential discontinuity using subprogram INVER 
to invert A matrix. 

•Computes C (flutter airforce) matrix. 
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CX:0N5,6 and CCONI 1,12 

•Computes terms in C-matrlx. 

ALIST: 

•Used in printing out various matrices. 

TEST: 

•Used in retrieval of Data vector for ACALC and CMATR. 

INVER: 

•Used In Inverting A matrix. 

Flow Diagram : A functional flow diagram of the computer program Is 

pictured In Figure A-l. The aerodynamic Influence matrix calculated by ACALC 
is inverted and used to compute generalized forces due to specified motions 
of the surfaces by CMATR and the appropriate CCONs. The use of generalized 
forces was found to offer considerable flexibility. For example, it Is 
possible to calculate the flutter airforce matrix or the lift slopes merely 
by specifying the modal displacements in the correct manner. 

DATA Vector ; The DATA vector contains the complete physical description 
of the airfoils and of the aerodynamic elements. Each block of data Is desig- 
nated by a "Type” number, with provisions for 15 types stored sequentially. 

A mixed Integer and floating point format Is used with the Integers coming 
first In each line. 

The first 60 words of the vector contain information required to identify 
the 15 types of data stored, and the 61st word represents the total number 
of lines. The complete vector contains the following: 

Identification: 

NADD( I ) = First address of Type I data (i = I,. ..,15). 

NFIX(i) = Number of integer words per line of Type I data. 

NFL(i-) = Number of floating words per line of Type i data. 

NL(i) = Number of lines of Type I data. 
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Fig. A-1 Flow Diagram 
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Type I ; Description of Airfoil Surfaces 
NSURL, NSURK = Surface number 

I I if surface is symmetric about z-ax Is 

-I if surface is symmetric about z-axis, but motion Is 
antisymmetric 

0 if there Is.no symmetry 
NSYMY = same as above, but relative to y-axis 

NGEOM = number of like surfaces on which aerodynamic forces are 
ca I cu I ated 

XO = x-origin of surface 
YO = y-origin of surface 
ZO = z-origin of surface 
COSG = cos g 
SING = sin g 

Type 2 : Point Coordinates 

L, K = Sending or receiving point number, I or k 
NSURL, NSURK = Surface number 
XL, XK = x-coordinate 
SL, SK = s-coordinate 

Type 5,7,11; Aerodynamic Element 

L = Sending point number, 1 
XFOR = Most forward side, XpQp 
XAFT = Most rearward side, X^pj 
SIN = Inboard edge, Sj^ 

SOUT = Outboard edge, SQy-j- 

Types 6, 8, 12 : Aerodynamic Element (Trailing edge) 

L = Sending point number, I 

NBOXES = Number of regions to be assumed in wake 
XFOR, XAFT, SIN, SOUT as above 
Element Types 5,6: Mid Point Constant Potential (MPCP) 

7,8: First Order, Constant Potential (FOCP) 

11,12: First Order, Zero Pressure (FOZP) 
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A Vector: In order to minimize the use of the core storage, the matrices 

In complex form are all stored In a one dimensional array Afl), Table A-l 
shows the composition of the A vector. The matrices stored. In All) are 
described as follows: 

A -Matrix 

The A-matrIx contains the comp lex, aerodynamic matrix A(k, £), and Is 
stored In complex form in the A vector with a blank column to facilitate 
matrix Inversion. The general term Is A(L) where L = (k - I) (NRP + 1) + £. 
The required total storage Is 2NRP(NRP +1). 

W~Matrix 

The W-matrix contains the modal values of the angle of attack, a(l, k), 
and deflection, h(l,k), at each point k due to mode I stored in the complex 
form a(I, k), h(i, k). W-matrIx Is stored In complex form In the A vector. 

The general term of W-matrIx Is ACNADDW + L) where L = (MODE - I) NRP + k 
and MODE takes the values I to (NMODEI + M0DE2). NMODEI Is the number of 
values of column Index I In the flutter airforce matrix C(I, j), and NMODE 2 
Is the number of values of column Index j in the C-matrIx. If NM0DE2 Is 
equal to zero, j belongs to the same set of modes as I, and C Is a square 
matrix as used In flutter calculations. This feature simplifies the calcu- 
lation of lift derivatives and distributions. The total storage required for 
W-matrIx Is 2NRP(M0DEI + NM0DE2). 

P-Matrix 

The P-matrIx contains the complex matrix of potential discontinuities 
P(j, k) and Is stored In complex form In the A vector. Its general terms 
Is A(NADDWF + J3) where NADDWF = ((NRP + I) + NMODEI + 2NM0DE2) NRP and J3 = 
NRP(MODE-I) + k. MODE takes the values I to NTEMP where NTEMP = NM0DE2 and If 
NM0DE2 = 0, NTEMP— NMODEI. Total storage required Is 2NRP( NTEMP). 

C-Matrlx 

The C-matrIx contains the complex flutter airforce matrix C(l,j) and Is 
stored In complex vector form in the matrix A. Its general term Is 
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Table A-1 Composition of A Vector 


Address 


A(I) 






NTOT 

Matrix 

Aerodynamic 

Matrix 

A(k,P) 

or 

its Inverse 

Modal 
Matrix 
w(i, k) 
for 

NMODE 1 

Modal 
Matrix 
w(i, k) 
for 

NMODE 2 

Same as 
w(i, k) 
for 

NMODE 2 
except 
Imag . part 
is multi- 
plied by 
Freq. ratio 
(to/v) 

Potential 
Discont. 
matrix 
P(jr k) 

Flutter 
Airforce 
matrix 
C(i, j) 



NADDW « (NRP + 1) NRP 
NADDW2 = NADDW + NRP X NMODEl 
NADDW3 = NADDW2 + NRP X NM0DE2 
NADDWF = NADDW3 + NRP X NM0DE2 
NADDC = NADDWF + NRP X NMODE2 
NTOT *= NADDC + NMODEl X NMODE2 


NJ 

04 


Appendix A 









Ill II 


■ III ■■mill ■■ III 


IIWI I III 


Appendix A 

A(NADDC + NTEMP(i-l) + j) where NADIX = NADDWF + NRP*NTEMP, and NTEMP, i and 
j take the same values as above. The total storage required is 2 NMODEI (NTEMP) . 

Symmetry : For each surface, the condition of symmetry is expressed by 

the two numbers NSYMZ and NSYMY, each having three values. In each case, 
a nonzero value indicates structural symmetry (the sign Indicates whether the 
downwash normal to the surface is positive or negative, respectively). Thus, 
there can be four conditions of structural symmetry corresponding to four 
values of a counter MSYM as shown in Figure A~2. For each value of MSYM, values 
are assigned to the four quantities KSYMZ, KSYMY, LSYMZ, LSYMY as shown. However, 
when the condition Indicated Is not met, the counter MSYM skips a value. 

To simplify computation, contributions to the A-matrlx are obtained by 
assuming the receiving points to be rotated into the four segments In turn. 

The sending points remain In the principal segment. This is accomplished by 
altering the equations developed In Sections II and III by multiplying cos g 
sin g, and AA by the following factors; 

cos g by KSYMZ 

sin g by KSYMY 

AA by LSYM = LSYMZ* LSYMY* KSYMZ 

Some examples of possible airfoil combinations together with the appro- 
priate values of NSYMZ and NSYMY are shown in Figure A-3. It should be noted 
that the symmetry conditions are supplied separately for each airfoil. 

Preparation of Input Data : The Input data are read by the main program 

FLUTT and by subprograms DMATR or DREAD, and WMAT or WREAD. The user can choose 
either OMATR or DREAD for input of Data vector and WMAT or WREAD for input of 
W matrix. For rectangular elements obtained by breaking up each surface Into 
chordwise strips of varying width and breaking up each strip into equal size 
rectangular elements, one can use DMATR by supplying the coordinates of each 
strip and the number of elements to be contained in each strip. The W matrix 
for surfaces divided into such rectangular elements can be easily generated 
with the use of WMAT by supplying the modes of motion for each strip. To 
divide a surface into rectangular elements of nonuniform size, the user should 
calculate all the required data for D vector and then read them In by . 

DREAD. The correspond! ng W matrix should also be calculated by the user and 
then read in by the subprogram WREAD. 
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MSYM = 3 

condition: NSYMZ =1,-1 

NSYMY =1,-1 

KSYMZ » -1 
KSYMY - -1 
LSYMZ - NSYMZ 
LSYMY -NSYMY 




MSYM = 2 

condition: NSYMZ =1,-1 

KSYMZ = -1 
KSYMY = 1 
LSYMZ = NSYMZ 
LSYMY = 1 


MSYM = k 

condition: NSYMY =1,-1 


KSYMZ = 1 
KSYMY = -1 
LSYMZ = 1 
LSYMY = NSYMY 



y 


point 



t 


MSYM = 1 

condition: none 

KSYMZ = 1 
KSYMY = 1 
LSYMZ = 1 
LSYMY = 1 


sending point 


Multiply cos g by KSYMZ; y , by KSYMZ 

ok 

sin g by KSYMY; Z^^ by KSYMY 

AA by LSYM = LSYMY*LSYMZ»KSYMZ 


Figure A-2. Symmetry Conditions 
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(a 


Note: 



) Wing with Dihedral 


NSYMZ = 0 
NSYMY = 1 


I 

I 


I 


\ 


z 


y 



(b) Symmetrical Vertical Tail 



NSYMZ = 

[ 

= 1 


NSYMY = 

^ -1 





z 


(c) Cruciform Tail (d) Biplane 

Figure A-3. Examples of Configurations 


Surface for which data is supplied 



Surface implied by symmetry conditions 
Direction of downwash 
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Before attempting to make use of the computer program, one should make the 
following selections: 

•Selection of subprograms to be used: 1) Input for D vector — DMATR or 

or DREAD; 2) Input for W matrix — WMAT or WREAD; 3) Calculation of A 
matrix — Select aerodynamic elements to be used. 

•Desired print out of various matrices. 

•Local coordinates for each surface. 

•Planforms for each surface. 

•Conditions of symmetry for geometry and motion of each surface. 

•Numbers of collocation points and chordwise strips. 

•Effective wake length. 

•Number of repeated runs using same D vector and W matrix for different 
Mach numbers and frequency ratios. 

•Mach numbers and frequency ratios. 

The user should also be aware of the inherent limitations in the program 
as follows (as set by the DIMENSION statements): 

'Number of collocation points must be < 72. 

•Number of chordwise strips must be < 14. 

•(Number or modes) X (Number of chordwise strips) must be < 84. 

•Number of surfaces must be <4. 

•Number of repeated runs must be < 10. 

**» Data Block F-0* »* 

Number of cards = I; Format (lx, 79Hbbb b) 

The above Is a sample of a title card. Anything placed on this card will 
be printed as output at the beginning of each run. 

*» *Data Block F~ l *** 

Number of cards = I; Format (20 13) 

This block specifies the print out of D, W, A, and P matrices. 

I : print, 0: no print 


I 
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[■*5 4~»6 7->9 I0 » I2 

( IC NPD NPW NPA NPP 

where IC = Card Identification number = I 

NPD = PrInt-out of D vector 
NPW = Pr I nt-out of W matrix 
NPA = PrInt-out of A matrix 
NPP = PrInt-out of P matrix 

*»* Data Block F-2* ** 

Number of cards = 1; Format (2013) 

This block shows the number of Integer words per line of Type I data. 
Types not used should be equal to zero. 

(I = I,. ..,15). 

1-^5 4-»6 37-»39 

f IC NFIX I (I = I, 2, ..., 15) 

where IC = 2 
NFIXI = 4 
NFIX2 * 2 
NFIX3 = 0 
NFIX4 = 0 

NFIX5, 7, II, = I for the type i used 
NFIX6, 8, 12=2 for the type I used 

*** Data Block F-3* ** 

Number of cards = I; Format (20I3‘) 

This block shows the number of floating words per line of Type I data. 
Types not used should be set to zero. 

I ->3 4->6 37-*»39 

( IC NFL I (I = I, 2, ..., 15) 

where IC = 3 
NFLI = 5 
NFL2 = 2 
NFL3, NFL4 = 0 

NFL5, 7, II, = 4 for the type I used 
NFL6, 8, 12=4 for the type I used 
128 
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*** Data Block F-4^ ** 

Number of cards = 1; Format (2013) 

This block shows the number of lines for each type of data. 

I ->5 4->6 37-»39 

(iC NLi (i = I, 2, 15) 

where IC = 4 

NLI = Number of surfaces 
NL2 = Number of collocation points 
NL3 = 0 
NL4 = 0 

NL5 to NLI 2: Equals the number of elements of the Type 5 to 12, 

respectively. 

*** Data Block 

Number of cards = I ; Format (2013) 

This block determines the selections of subroutines for D vector 
and W matrix and the number of repeated runs to be made with the 
same D vector and W matrix for different Mach numbers and frequency 
ratios. 

I ->-5 4-»^ 7->9 i0->-l2 

lie NDMAT. NWMAT NAMAT 
where IC = 5 

NDMAT =1 if subprogram CM^ATR is used and 2 if subprogram 
DREAD is used. 

NWMAT = I if subprogram WMAT is used and 2 if subprogram 
WREAD is used. 

NAMAT = Number of runs to be made using the same D vector 
and W matrix for different Mach numbers and 
frequency ratios. 
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** *Data Block 

Number of cards = I; Format (2013) 

This block specifies the following items; 

4-^ 7->9 1012 13^13 : 16^18 

fiC NMODEi NM0DE2 NRP MTOT NTE 

where IC = 6 

NMODEI = Number of values of row Index I In C matrix 
NM0DE2 = Number of values of column Index j In C matrix 
NRP = Number of reference points (col location points) 

MTOT - Number of chordwise strips 
NTE = Number of trailing edge elements 


Data Block DM-7 to Data Block DM-22 are used only when NDMAT of Block F-5 

is I . _ 

** *Data Block DM-7 *^* 

Number of cards = 1; Format (2013) 

K3 4-^6 

fiC (NSUR(I), 1=1, NL) 

where IC = 7 

(NSUR(I), 1=1, NL) = Surface number for each surface 
NL - Number of surfaces 
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** *Data Block DM=8* ** 

Number of cards = I; Format (2013) 

1^5 4-> 6 

f IC (NSYMZ(I), I = I, NL) 

where IC = 8 

(NSYMZ(I), 1=1, NL) =1 if surface is symmetric about z-axis 

= -I if surface is symmetric about z-axis, but 
motion is antisymmetric 

= 0 if there is no symmetry 

*^* Data Block DM-9* ** 

Number of cards = I ; Format (2013) 

I ->3 4->- 6 

pc (NSYMY(I), 1=1, NL) 

where IC = 9 

(NSYMY(I), 1=1, NL) = same as in Block DM-8, but relative to y-axis 

*** Data Block DM- 10*^ ^ 

Number of cards = I; Format (2013) 

I ■>•3 4"^ 6 

(jC (NGEOM(I), 1=1, NL 

where IC = 10 

(NGEOM(I), 1=1, NL) = Number of like surfaces on which aerodynamic 

forces are calculated 

*^^ Data Block DM- 1 1^ *^ 

Number of cards = I; Format (IX, 12, 4X, 13, 6FI0.5) 

l-v3 8->l0 I K 20 

[ IC NC ~ (X0(I), 1=1 , NL) 

where IC = I I ; NC = Card number in each block (i.e., for 1st card NC = 1 ) 
(X0(I), I =' I , NL) = x-origin of surface I 


131 



Appendix A 


- ! o ck DM- 1 2* ** 

of cards = I; Format (IX, 12, 4X, 13, 6FI0.5) 

I ^>3 8^10 I K 20 

NC (Y0(I), 1=1, NL) 

IC = 12; NC = Card number in each block 
.'MI), I = I, NL) = y-origin of surface I 

''-H ock DM- 1 5* ^^ 

or of cards = I; Format (IX, 12, 4X, 13, 6FI0.5) 

!-> 3 8-^10 I I 20 

f rC NC (Z0(I), 1=1, NL) 

IC = 13; NC = Card number in the block 

~j(I), 1=1, NL) = z-origin of surface I 

-0 Bl ock DM-14** * 

of cards = I; Format (IX, 12, 4X, 13, 6FI0.5) 

- 8-^10 I 1^ 2Q 

NC (COSG(I), 1=1, NL) 

0 .) =14; NC = Card number 
J(I), 1=1, NL) = cos g of the surface I 

■ ok DM- 1 5* ^^ 

of cards = I; Format (IX, 12, 4X, 13, 6FI0.5) 

H MO IK 20 

’in (SING(I), 1=1, NL) 

i IC = 15; NC = card number of the block 

^’SINGd), 1=1, NL) = sin g of the surface I 
g angle between local normal and z-axis, positive clockwise 

ck DM- 1 6* ** 

. of cards = I ; Format (2013) 

4-^6 


(NSURF (I), 1=1, MTOT 


Appendix A 


where IC = 16 

(NSURF (I), I*=l, MTOT) = Surface number of strip I 

** *Data Block DM- 1 7* ** 

Number of cards = I; Format (2013) 

|«»5 4~*~6 

^^"iC (NQT(I), 1=1, MTOT) 

where IC = 17 

(NQT(I), I - I, MTOT) == Number of elements In strip I 

** *Data Block Dh4-I8* ** 

Number of cards = I; Format (2013) 

K3 4 -» 6 

(NB(I), 1=1, MTOT) 

where IC = 18 

(NB(I), 1=1, MTOT) = Number of elements to be assumed In the wake behind 
strip I 

** *Data Block DM- 1 9* ** 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

K3 &»I0 II ->-20 

(c NC (SF0R(I),I = I, MTOT) 

where IC = 19; NC - Card number In the block 

(SFOR(I), I = I, MTOT) “ Leading edge of strip I 

** *Data Block DM-20^ ** 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

1-»3 8^10 II 20 

C NC (SAFTtt), 1=1., MTOT) 

where IC = 20; NC = Card number In the block 

(SAFT(I),I = I, MTOT) = Most rearward edge of strip I 

»* *Data Block DM-2}* ** 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

l->-3 8-»l0 llH-20 

NC (SIN(I), 1=1 , MTOT) 
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where, IC = 21; NC = Card number In the block 
(SIN(I), I='i , MTOT) - Inboard edge of strip I 

^** Data Block DM-22^ ** 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6F10.5) 

I ~>-5 1 0 I l~^2Q , 

[^~IC NC (SOUT(I), 1=1, MTOT) 

where IC = 22; NC = Card number in the block 

(SOUT(I), 1=1, MTOT) = Outboard edge of strip I 


Data Block DR-7 to Data Block DR- 10 are .used only when NDMAT of Block F-5 
is 2. 

^** Data Block DR-7 and DR-8* *^ 

These blocks supply type I data. Integer words are read on cards DR-7 
and floating words on cards DR-8. Repeat DR-7 and DR-8 for each 
surface. 

DR-7 ; Number of cards = Number of surfaces; Format (2013) 
l->3 4->6 7->9 ICh-12 15^15 I6->I8 
[iC IT NS UR NSYMZ NSYMY NGEOM 

where IC = 7 
IT = I 

NSijR = surface number 
NSYMZ = See Data Block DM-8 
NSYMY = See Data Block DM-9 
NGEOM = See Data Block DM- 10 

DR-8: Number of cards = Number of surfaces; Format ( IX,I2,7X,6F10.$) 

l->5 4->6 II -» 20 21 30 51 40 41 50 51 -» 60 

icT XORIG YORIG ZORIG COSG SING 

where IC = 8 

IT = I 

XORIG = x-origin of surface 

YORIG = y-origin of surface 

- z-origin of surface 

COSG = cos g 
SING = sin g 
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** *Data Block DR-9 and DR- 10 *** 

These blocks supply type 2 data. Integer words are read on cards DR-9 
and floating words on cards DR- 10. Repeat DR-9 and DR-10 for each 
collocation point. 

DR-9 : Number of cards = Number of collocation points; 

Format (2013) 

K3 4-»6 7->9 10 12 

IT PT NSUR 

where IC = 9 
IT = 2 

PT = Collocation point 
NSUR = Surface number 

DR-10: Number of cards = Number of collocation points; 

Format (IX, 12, 7X, 6FI0.5) 

l->-3 4->^ 11^20 2l->-30 

(iC IT X S 

where IC = 10 
IT = 2 

X = x-coordlnate of collocation point 
S = s-coord i nate of coll ocat Ion po I nt 

*** Data Block DR“I I and DR- 1 2* *^ 

These blocks supply the data for the on-the-wing aerodynamic elements. 
Integer words are read on cards DR- 1 I and floating words on cards 
DR- 12* Repeat DR- 1 I and DR- 1 2 for each elementi 

DR- 1 1 : Number of cards = Number of elements 
Format (2013) 

I •♦■3 

(jC IT L 
where IC - 1 1 

IT = 5, 7, or 1 1 depending on the selection of aerodynamic elements 
L = Sending point number 

DR- I 2 : Number of cards = as for DR-II; Format (IX, 12, 7X, 6FI0.5) 

l-»-3 4->^ I l->20 ’2l-»30 3l->40 4l-»50 
IT XFOR XAFT SIN SOUT 
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where IC = 12 

IT = As i n DR- 1 1 

XFOR = Most forward side of an element CX-Coord.) 

XAFT = Most rearward side of an element 
SIN = Inboard edge of an element (S-Coord.) 

SOUT = Outboard edge of an element 

** *Data Block DR- 1 3 and DR- 1 4* ** 

These blocks supply the data for trailing edge elements. 

Integer words are read on cards DR- 1 3 and floating words on cards 
DR- 1 4. Repeat OR- 1 3 and DR- 1 4 for each element. 

DR- 1 3 : Number of cards = Number of trailing edge elements; 

Format (2013) 

K5 4-»6 7-»9 10^12 

^IC IT L NBOXES 

where IC = 1 3 

IT = Type of data 
L = Sending point number 

NBOXES = Number of regions to be assumed In the wake 

DR- 1 4 : Number of cards = as for DR-13; Format (IX, 12, 7X, 6FI0.5) 

l->-3 4-»^ II 20 21 30 31-^ 40 41 50 

IT XFOR XAFT sTn SOUT 

where IC = 14 

IT = As In DR- 1 3 

XFOR = Most forward side of a trailing edge element 
XAFT = Most rearward side of wake region 
SIN = Inboard edge 
SOUT = Outboard edge 


Data Block WM-I to Data Block WM-6 are used only when NWMAT OF Block F-5 
is I . 

*^ *Data Block WM-I* ** 

Number of cards = NMOMT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

l-»3 8-*-IO I I ->• 20 

j*IC NC (H(I), 1 = 7 , NMOMT) 
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where IC = 23; NC = Card number in the block 

I = I, NMOMT) = Deflection of strip K due to mode J 
I = ((J-l )*MTOT t K) 

J = Mode number 
K = Strip number 

NMOMT = CNMODE I + NMODE 2)*MT0T;Clf NM0DE2 = 0, NMOMT = NMODE I *2*MT0T] 
** *Data Block WM-2* ** 

Number of cards = NMOMT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

K3 8->IO ll-)-20 

NC (ALPH(I), 1=1, NMOMT) 

where IC = 24; NC = Card number in the block 

(ALPH(I), 1=1, NMOMT) = Rotation about XROT for strip K due to mode J 
See Block WM-I for other definitions. 

*** Data Block WM-3* *^ 

Number of cards = NMOMT/6; Format (IX, 12, 4X,. 13, 6FI0.5) 

l-»3 8->IO I I 20 

(^C Nc"^ (XROT(I), 1=1, NMOMT) 

where IC = 25; NC = Card number in the block 

(XROT(I), 1=1, NMOMT) = Axis of rotation for strip K due to mode J 
See Block WM-I for other definitions. 

** *Data Block WM-4^ ** 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6F10.5) 

l->-3 &»I0 I l-> 20 

(k NC (SIN(K), K=l, MTOT) 

where IC = 26; NC = Card number In the block 
(SIN(K), K=l, MTOT) = Inboard edge of strip k 

** *Data Biock-WM-5^ *»^ 

Number of cards = MTOT/6; Format (IX, 12, 4X, 13, 6FI0.5) 

K3 \\^20 

I^C NC (S0UT(K),K=l , MTOT) 
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where XC = 27; NC = Card number In the block 

CSOUTCk), K=l, MTOT) = Outboard edge of strip K 

*** Data Block WM-6^ ** 

Number of cards = I; Format (2013) 

l-»5 4 ->■ 6 

flC (NSUR(K) (K=l, MTOT) 

where IC = 28; NC = Card number In the block 

(NSUR(K), K=l, MTOT) = Surface number of strip K 


Data Block WR-I Is used only when NWMAT OF Block F~5 Is 2. 

*** Data Block- WR-1* ** 

Number of cards = NRP*NM0DE*2/6; FormatdX, 12, 4X, 13, 6FI0.5) 

J-^3 &^\0 I I. -^20 

NC CA(I), 1=1, NTOTW) 

where IC = 23; NC = Card number in the block 
(A(I), I- I, NTOTW) = Ith element of W matrix 
NTOTW = NRP*NM0DE 

** *Data Block F-7* ** 

Number of cards = NAMAT/6; FormatdX, 12, 4X, 13, 6FI0.5) 

l-)-3 8->IO II ->>20 

NC (RM(I), 1=1, NAMAT) 

where IC = 29; NC = Card number In the block 

(f^(I), 1=1, NAMAT)=Mach number for Ith repeated run 

** *Data Block F-8* »* 

Number of cards = NAMAT/6; Format(IX, 12, 4X, 13, 6FI0.5) 

|->-3 8->-IQ I I 20 

|IC NC (FRA(I), 1=1, NAMAT) 

where IC = 30; NC = card number in the block 

(FRA(I), 1=1, NAMAT) = Frequency ratio for Ith repeated run 
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Example Problems. - 

Example I: Input data for the example In Figure A-4 Is read by 

using subprograms DREAD and WREAD for rectangular elements of nonuniform size. 
In calculating A and C matrices, AC0N5 and AC0N6 and CC0N5 and CC0N6 were used. 
The outp shows the listing of D vector, W, A, P, and C matrices. 

Definitions of Modes and Aerodynamic Coefficients: The oscillating modes 

considered - 

J (X) =-( * , n = 1 

" 'X-XROT , n = 2 

The aerodynamic coefficients considered - 


C. . = // J (X)AP^(X)dXdy 
i,j m - n - ^ 



Mode 

m, n 

a 

h 

^RP 

^ LU 

XI Q 
— O 

E S 

~~3 Z 

— » CM 

XI UJ 
Q 

c o 

-5 2 

z 

1 = 1 

m=l 

0 

2.0 

K5 

1=2 

m=2 

4 

Sc 

|<X-.5) 

K5 

■ 

n=2 

1.0 

X-.5 

K5 


S = Half Wing Area = 1.0 
c = Root Chord = 1.0 
AP^ = Pressure jump corresponding to 

Example 2: The example Illustrated In Figure A-5 shows how the Input 

data for- a Swept Wing with Partial -Span Flaps can be read In by using sub- 
programs DMATR and WMAT. The A and C matrices were calculated by using AGON 1 1 
and AC0NI2 and CCONII, CC0NI2, respectively. Portions of various output 
listings are shown for reference. 
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Fig. A-4 Rectangular Wing of Aspect Ratio Two, 
Pitching about Midchord 



Append tx A 


Input for Example J 


Input 


tXAr->HLL‘is 

T<nrrrTimr 

USl^C DK^AU AfjU WRLAU 



1 

1 

1 

1 

1 








H 

2 

□ 


QIIQJQ] 

-0 -0 -0 - 


-0 -0 



s 

5 

2 

0 

0 4 4 

PBIHI 

-0 -0 -0 - 

0 

-0 -0 



k 

1 

5 

U 

U 3 2 

-0 -u -u 

-0 -0 -0 - 


-0 -0 



9 

2 

2 

1 








6 

2 

1 

9 

0 







7 

1 

1 

I 

0 1 







6. 


1 


u.ooooo 

o.oooou 

0.00000 


l.OOOQO 

0.00000 


9 

2 

■ 1 

1 








10 


2 


•37500 

.10000 







2 

2 

1 








10 


2 


.37500 

.50000 






9 

2 

3 

1 








10 


2 


•37500 

.90000 







2 

4 

1 








10 


2 


•875u0 

•25000 







2 

5 

1 








10 


2 


•B7500 

•75000 






11 

5 

1 









12 


5 


•ipsoo 

.62500 

0.00000 


•20000 



11 

5 

2 









12 


5 


•1?500 

.62500 

•20000 


•eoooo 



11 

5 

3 









12 


5 


•12500 

.62500 

.60000 


1.00000 



13 

o 

4 

1 








IH 


6 


•62500 

1.12500 

0.00000 


• 50 ono 



iS~ 

6 

5 

1 








It 


6 


•62500 

1.1250U 

•50000 


1.00000 



23 


1 


o.onooo 

2.0UO00 

0.00000 


2.00000 

0.00000 

2.00000 

23 


2 


0 • 0 p 0 C u 

2 , «) fi t) U U 

o.unooo 


2.00000 

4.00000 

-.50000 

23 


3 


4 .OoUtU 

-.50000 

4.00000 


-.50000 

4.00U00 

1.50000 

23 


4 


4 .OoCliU 

1.50U0G 

l.GOOOO 


-.12500 

1.00000 

-.12500 

23 


5 


1 • 0 r 0 0 u 

-.12500 

l.ODC'OO 


.37500 

1. 00000 

.37500 

"29' 

— 

— 1“ 

: 








30 


1 


o.onocu 

-*f— — V """ ■ 
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.Output for Example I 

proSrah for unsteady air forces, using 

oownwash-velocity potential method 

BY J.K. HAVILAND AND Y.S. YOO « UNIV OF VA 
PART 1 • READ IN W AND DATA MATRICES 


EXAMPLE! I RECT*. WING USING DREAD- AND WREAO 

11111 

2- — — 2- — 0 — fl - --l — 2 -*0_*0 .“0 .-O.-rO -•0 - *0. 

3 S 2 0 0 4 4 -0 -0 -0 “0 -0 -0 -0 -0 -0 
— 5- 0 -- 0- 3 — 2--0 --0--0 --0 --0 --0- -0 -0— f O- 
5 2 2 1 


— 6 - 2 -15 - 0 -0 

DATA matrix ELEMENTS 


. 62- 

-.4.. 

5 


71. 

_ 2 

-.2.! 

5. . 


0 

0 

0 - 

91 

0 

0 

0 

91 

1 

4 

3 

106 

2 

4 

2 

lie 

-0_ 

— -0_ 

_-0. 

- 116 

•• Q 

_-0 

.-0 

- 118 

. -0— 

.-0„ 

_-0 

116 

-0 

-0 

-0 

116 

-0 

-0 

-0 

116 

-0 

-0 

-0 

116— 

• 0.. 

— -0- 

•0-^ 

— 116. 

- -0. 

__-Q. 

_-0 

-118- 

„-Q— 

-0_ 

-0 


118 

1 — 1 1 0._1 


8 

9. 

2 

1 

X 1- 

0.00000 

0.00000 

0.00000 

1.00000 

0,00000 


10 

A 


2 

,37500 

• 10000 





10 

9_ 


2 

^ 1 

•37500 

•50000 





10 

9. 


2 

tt /f 

•37500 

•90000 





10 

9. 

2 

9 ^ 

•87500 

• 25000 





10 

11 

•s 

2 

1 

•67500 

•75000 





12 

-11 

5 

5 

2 

•12500 

• 62500 

0.00000 

•20000 



12 

11 

12 

-13- 

14 

13 

14 
23 
23 
23 
23 
.23_ 

5 

5 

.3 -. 

•12500 

.62500 

•20000 

•80000 



6 

5 

4 1- 

•12500 

0 

0 

tr> 

CM 

Ml 

• 

•80000 

1.00000 



6 

6 

5—1 

•62500 

1.12500 

0.00000 

•50000 




6 

-1- - 

•62500 

0.00000 

1.12500 
2.00000 . . 

•50000 

0.00000 

1.00000 

2.00000 

0.00000 

2,00000. 


2 

3 

4 

_5 

0.00000 

4.00000 

4.00000 
_1.00000_ 

2.00000 
• •50000 
1.50000 
-.12500 

0.00000 

4.00000 

1.00000 
1.00000_. 

2.00000 
-.50000 
-.12500 
^37500 

4.00000 

4.00000 

1.00000 
__l.Q0000_ 

-.50000 
1,50000 
-.12500 
__ .3750(1 
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Cont’d Output for Example I 


WMATKIX tLEKLliTS 


ROW COL i^LAL IMAGINARY 

1 1 0. 2.OU0OOOOUE + O0 

120. 2 , UOOOOOOUE+00 

1 3 .. C. . 2.0000000UE.+ 00 . 

1 H 0. 2.UOOOOOOUE+00 

1 5 1 . 2 . 0 0 0 0.0 0 0 U E + 0 0_ 

2 1 H .OOOUOOuOE + 00 -5, OOOOOOOOE-Ol 

- .2 2 .4 . 000000U0E+00_... -5. OOOOOOOUE-Ol 

2 3 4, OOOUOOUOE+UO -b. 00C0000UE~0l 

. 2. .4 ^.QOOUUOuO.LtUQ l.bOOOOOOUE + 00 

2 5 4 , OOOOOOUOe + 00 1 . 5000000UE + 00 

3 1 l«QO.O.uGO.UOEtO.O . -1.2b 0.0 OOOUE-Oi 

3 2 1 .000UUGU0E4-00 -1 . 2500000UE'Ul 

3 3 1 .OOOUOOOOE + 00 _ - 1 . 250 0 0 0 0 UE" 0 1 

3 4 1 . OriOuOOUOL + 00 3.7500000UE-01 

.3_.. .5 1 .OOUUOOUOE + 00 3.7500000UE-01 

29 1 L.OqUUU 

_3P 1 O.COOUO 


■PA1TY'"2"‘ -njTKlVE 7\ERonYuAMlC MAlTTlX 

A M AT R I X ■■ ■ 1‘ f/KH= ■ "b^PiACH NO^OTOT) 0 C'O F”K^:r“07W0 0 0' 

ROW COL plal imaginary 

1 1 in 

1 2 -1 .2nGj.i /b3L + CG 0. 

I — 3 — o ~ 

1 4 -3.2aRo‘+234l-01 0. 

1 b"" -r.5r570 5^3‘9L-'CJ2‘ ^0“^ 

1 G 0 • 0 . 

2 1 - 1 . b ? bb u 04 L - O'l ir> 

2 2 1 .b?9oL.fJ7bE + U0 0, 

2 T -r.-T^lH^bOI-in IT. 

2 '♦ -1 .954^7 ?1 95L-Ui 0. 

2 5 - I . T 3 l 9 TIX 3 L ‘^ in : — 

2 b U . 0 • 

3 ‘ 1 “.2.3i479Bi?L^C;2 07 

3 2 -1 .047G3G37E + 00 0. 

3 3 ~3'.'4 2 b 9 4 S'b 0 1 +“^0 ' “ 

3 4 -5.307ti9BlbE-02 0. 

3 ■ ' 5 ' 7 2.0A ib7‘+ Vll> Cr^ U7 ” 

36^0. 0 . 

■ 4 1 7lV2^19bb3E^^1Cri 07 

4 2 -2.37l‘ib3y4E-01 0. 

'““4 3“ .93nb3241L-02 7“^' “ ” 

4 4 ' 1.57137243L+00 0. 


I 


H3 


Cont*d Output for Example I 


\ 


u 

5 

-«;*5lllti69dUL-0l 

U. 


4 

6 

Lr • 

0. 



1 

-k- .37-i^7l3RL.02 

0. 


5 

? 

.0Afo3271L“0l 

a. 


S) 

3 

.L,Pi9du7bkL;-Uk 

0. 



4 

-k,b;utia'3dOL-oi 

c. 


b 

b 

1 ,7‘^‘»4ulii'5L + 0U 

u. 


5 

6 

C'. 

0. 


PAHT 

3- 

DLKIv/L CVAIKIX 




list; 

irjo up 

PUTFf.UAL UiSCON { lUUiTlt.S 

ROt^i 

COL 

Rl AL 

IMAGIlvAKY 

1 

1 

1 .OllioUlOl L + LC 

0. 

1 

2 


0 • 

1 

3 

t ,f,;^72‘:)7i2L“Ul 

0 . 

1 

4 

1.0] 13o540E+0U 

0. 

1 

5 

P.S7t)ldl 7P[1-01 

L' • 


C*'VA7 jriX ' F O'K li= 1 


MACH 

li 

• 

o 

{'..UOODU FKA= 0.00000 

ROw 

COL 

AL. IMAGIfvAKY 

1 

1 

1.06oU03b7L+00 0. 

2 

1 

— 1.3ob3/0cibL + 00 0. 


144 




Appendix A 



Control Surface Chord = .712 x half local chord 

Wing Aspect Ratio = 2 

Wing Taper Ratio = .2378 

Leading Edge Sweep Angie = 60 deg. 


Fig. A-5 Chordwise Strip Arrangement for a Swept V/ing 
Partial -Span Flaps 




J„(X) J (X) 

n — m — 

NM0DE2 NMODEl 
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Definitions of Modes and Aerodynamic Coefficients 


The oscillating modes considered; 


1 

For Wing 

9 

n"l 

X-Xjn 

For VJing 

9 

n=2 

1 

For Control 

Surface , 

n=3 

IX-Xc 

For Control 

Surface, 



The aerodynamic coefficients considered; 

AP„(X)axdy 

s 

where APjj = pressure jump corresponding to 
s = half span 
S = half wing area 


Mode I m, n 


m=l 


j=l n=2 


2 n=3 


1 -^ 


ALPH(I) 

Strips 

0 

l->6 

2 

S 

1-^-6 


l->6 


0 

1.5153 

1.7201 


1 -^ .732 l-s-6 .808 


0 1 ^ 

1.0 5-^6 


Strips 


j=3 n=U 


0 


0 

1.0 


l->-U 

5 ^ 


0 

1.5153 

1.7201 


1-i-li 

5 

6 



























































Appendix A 


InjDjjt for Example 2 


EXAMPLE PS <j]^EPT WINn KTTH CQfJTROL SURFACE, NRP=l({, 


1 

p. 

1 

4 

10 1 
2 -0 -0 -fl -0 

-0 -0 -0 

-012- 

n -n -0 



3 

5 

2 -0 -n -0 -0 

amoBa 

-0 4 4 - 




4 

2 

14 -0 -n -n -0 

-0 -9 -n 

-0 10 4 - 

0 -0 -0 



5 

1 

1 5 






6 

4 

3 14 6 4 






7 

1 

2 






A 

1 

1 • 






9 

0 

0 






10 

1 

!• 






11 


1 "OoOnooo 

-0.00000 





12 


1 -0,00000 

-0,00000 





13 


1 -o.oooon 

-0.00000 





14 


1 i.oonoo 

1.00000 





15 


1 -0,00000 

-0,00000 





lA 

1 

1112 2' 






17 

3 

2 2 3 2 2 






la 

16 

0 0 16 11 11 






19 


1 .34640 

.86600 

1.21240 

1.45880 

1.51540 

1.72010 

20 


1 1,71360 

1.51540 

1.72010 

2.05490 

1,85990 

1.95750 

21 


1 0.00000 

.40000 

.62000 

.84000 

.40000 

•62000 

22 


1 ,40000 

,62000 

.84000 

I.oonoo 

.62000 

• 84000 

23 


1 2,00000 

2.00000 

2.00000 

2.nonno 

2.00000 

2.00000 

23 


2 -0,00000 

-0,00000 

-0.00000 

-n.noooo 

-o.ooono 

-0,00000 

23 


3 -U.UllOuO 

-O.uOOuO 

“0,00000 

— U . j'j M fl U 0 

2 • u u u li u 

2 . OCOuO 

23 


4 -0.00000 

-0.00000 

-0.00000 

-o.nnooo 

-0.00000 

-0.00000 

23 


5 -0.00000 

-0.00000 

-o.ooono 

-o.nnono 

-o.ooono 

-0.00000 

23 


6 -0,00000 

-0.00000 

-0.00000 

-o.nnooo 

1,00000 

1,00000 

23 


7 -0.00000 

-0.00000 

-0.00000 

-o.nnooo 

-o.ooono 

-0,00000 

24 


1 -0.00000 

-QoOOOOO 

-0,00000 

-o.nnooo 

-o.ooono 

-0,00000 

24 


2 2.00000 

2.00000 

2.00000 

p.nnnoo 

2.00000 

2.00000 

24 


3 -O.OOQOO 

-0.00000 

“0.00000 

-n.oonoo 

-0.00000 

-0.00000 

24 


4 -0,00000 

-0.00000 

-0.00000 

-o.nnooo 

2.00000 

2.00000 

24 


5 o73200 

.73200 

.73200 

.732C0 

.73200 

.73200 

24 


6 -0,00000 

-o.ooono 

-0,00000 

-n.nonoo 

-o.ooono 

-0.00000 

24 


7 -0,00000 

-o.ooono 

-o.ooono 

-o.nnooo 

1,00000 

1.00000 

25 


1 -0.00000 

-o.oonoo 

-0,00000 

-o.nnooo 

-o.ooono 

-0.00000 

25 


2 .aoaoo 

080800 

.80800 

.80800 

,80800 

,80800 

25 


3 -OcOOboo 

-o.ooono 

-0.00000 

-o.oonoo 

-o.ooono 

-0.00000 

25 


4 -OoOOOOO 

-0.00000 

-0.00000 

-o.onooo 

1.51530 

1.72010 

25 


5 08O8OO 

csosno 

080800 

,80800 

.80800 

• 80800 

25 


6 -OqQOOOO 

-o.ooono 

-0.00000 

-o.nnnoo 

-o.ooono 

-0.00000 

25 


7 -0,00000 

-0.00000 

-0,00000 

-o.onooo 

1.51530 

1.72010 

2A 


1 0,00000 

•40000 

o620n0 

.84000 

•40000 

•62009 

27 


1 .40000 

,62000 

,84000 

1.00000 

•62000 

.64000 

28 

1 

1112 2 






29 


t 0,00000 

oSonoo 

,70000 

0 Annoo 

,90000 


30 


1 loOOOQO 

I.oonoo 

1.00000 

1 oOnnoo 

1.00000 
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Appendix A 


Ou+pu+. for Examp I e 2 

PROGRAM FOR UNSTEADY AIR FORCES USING 
OOHNWASH-VELOCITY POTENTIAL METHOD 

BY J.K. HAVILANO AND Y>S> YOO - UNIV OF VA 


PART 1 • READ IN H AND DATA MATRICES 
DATA MATRIX ELEMENTS 


62 

136 - 

136 - 

136 - 

210 - 
210 


1110 1 

•0*00000000 *0*00000000 *0*00000000 l.O'QOOOObO *0*000000 

1210 1 . 

1 *0*0 do on on 0 ' - o • ooodoooo odd o o do o i . odd do o oo -b.onDordo 


2 11 
2 *68820000 
2 2 1 

2 1*1439 333 3 

2 1*10952500 

2 4 1 

2 1*43422500 

2 5 1 

2 1*40278750 

2 6 1 

2 1*65663750 

2 7 1 

2 1*68282500 

2 8 1 

2 1*84819167 

2 9 2 

2 1*64458750 

2 10 2 

2 1*80912500 

2 11 1 

2 1.59966667 


0 • Ododoood odd o o do o i . d d o d d d oo -b.odOono 

•20000000 

.20000000 ^ 

*51000000 

- ^ ) i ■ 

. 51000000 
•73000000 
•73000000 
•92000000 

• 92000000 

•51000000 

•73000000 

*20000000 


2 

12 1 


2 

2*01355833 

•92000000 

2 

13 2 


2 

1,98908750 

*51000000 

2 

14 2 


2 

2*04652500 

•73000000 

11 

1 


11 

•46033333 

•91606667 

11 

2 


11 

*91606667 

1.37180000 


•40000000 

•40000000 
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Appendix A 

Cont*d Output for Example 2 


— MMATRliTELEMENT^ 


ROM 

COL 

REAL 

IMAGINARY 

1 

1 

-0. 

Z.OOOOOOOOE^OO 

1 

2 

-0. 

2«OOOOOOOOE^OO 

1 

3 

-0. 

2.00000000E'l>00 

1 

4 

-0. 

2.00000000E400 

1 

5 

-0. 

2.00000000E40Q 

1 

6 

-0. 

2,00000000E^00 

1 

7 

-O, 

2*00000000E^00 

1 

6 

-0. 

2.00000000E4>QQ 

1 

9 

-0. 

2.00000000E400 

1 

10 

-0 • 

2.00000000E4-00 

1 

11 

-0, 

2.00000000E4-00 

1 

12 

-0. 

2. OOOQOOOOE-fOO 

1 

13 

-0. 

2«OOOOOOOOE<fOO 

1 

14 

-0* 

2«0000OO00E-t-00 

2T 

1 

2.00 0DirODirE+0^ 

-2.396ddT^00E*0l 

2 

2 

2. OOOOOOOOE’t'OO 

6.71866667E-01 

2 

3 

2.0000000QE<t>QO 

6* 03050000E-01 

2 

4 

2.00000000E^00 

l«25245000Ei>00 

2 

5 

2.00000000E4>00 

1«16957900E^00 

2 

6 

2. 00000000E4-00 

1.69727500E4-Oa 

2 

7 

2 . dodoodroE^do 

1.74965dddE4>00 

2 

6 

2«OOOOOOOOE4>00 

2.08038333E'«’00 

2 

9 

2*O0000O0OE'('00 

l*67317500E<l-00 

2 

10 

2.00000000E40Q 

2«00229000E<»>00 

2 

11 

2.0000000QE4-00 

1.98333333E'f00 

2 

12 

2«00000Q00E*»>00 

2.41111667E<fOO 

2 

13 ’ 

2*OOOOOOddEtOd 

2*^62l75d0£<f00 

2 

14 

2.00000000E4-00 

2.4770S000E-f00 

3 

1 

-0. 

0. 

3 

2 


0. 

3 

3 

-0. 

0. 

3 

4 

-0. 

0, 

3 

5 

-0 • 

0« 

3 

6 

-0, 

0. 

3 

7 

-0. 

0. 

3 

6 

-0. 

0. 

3 

9 

-0. 

2.00000000E4-00 

3 

10 

-0* 

2. 00000000E4-00 

3 

“11 

-0. 

0« 

3 

12 

-0. 

0. 

3 

13 

-0. 

E.OOOOOOOOE-t-OO 

3 

14 

-0* 

2.00000000E4-00 

k 

1 

-0. 

0. 

4 

2 

-0. 

0* 

4 

3 

-0. 

0« 

4 

4 

-0. 

0. 

4 

5 

-0. 

0. 

4 

6 

-0. 

0. 
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Appendix A 


Cont*d for Example 2 


PART 

2 - 

DERIVE AERODYNAMIC 

MATRIX 

A MATRIX 

1 NRP= 14 MACH NO= 

•76000 FRAs 1.66660 

ROW 

COL 

REAL 

IMAGINARY 

1 

1 

1.47645344E+00 

7.40374286E-03 

1 

2 

-2.30407783E-01 

9.50228066E-02 

1 

3 

-4.92222006E-02 

2.38331129E-02 

1 

4 

-8.88720944E-03 

l,29734819E-02 

1 

5 

-5«91692197“E^3 

9.i956954iE-d3 

1 

6 

-7.72368777E-04 

5.73546826E-03 

1 

7 

-2.04522320E-05 

2.46429225E-03 

1 

6 

5,61022988E-04 

1.81960463E-03 

1 

9 

-9.91527271E-04 

4.87913843E-03 

1 

10 

3.61765554E-04 

2.18064490E-03 

1 

li 

-1,24052829E^04 

~3.34679343E-02 

1 

12 

3.22993542E-03 

3.31610288E-04 

1 

13 

4.02740113E-03 

4.02313751E-03 

1 

14 

5«36009830E*a3 

7,69676263E-04 

1 

15 

0. 

0. 

2 

1 

-2,04411530E-01 

-4.33756905E-02 

2 

2 

l,47645344E+00 

7.40374286E-03 

2 

3 

-2.96346732E-01 

2. 10531821E-03 

2 

4 

-1.01923316E-01 

2.55039169E-02 

2 

5 

-3.68973917E-02 

1.24130106E-02 

2 

6 

-1.34087399E-02 

l,04745443E-02 

2 

7 

-4.03207363E*03 

4. 14043549E-03 

2 

8 

-2.d29i2877£-03 

3.61443247E-03 

2 

9 

-1.51234926E-02 

1.05108725E-02 

2 

to 

-3.33496549E-03 

4.45918174E-03 

2 

11 

-2-15677709E-01 

1.25451330E-01 

2 

12 

3.26458742E-05 

5,51131065E-03 

2 

13 

-3,70683910E-05 

1.66086156E-02 

2 

14 

4.38441613E-03 

1.09612281E-02 

2 

15 

0. 

0. 

3 

1 

-1.46780441E-01 

-2. 72107623E-02 

3 

2 

-1.06648139E+00 

6.61679208E-03 

3 

3 

3.20433063E>0Q 

2.25398425E-03 

3 

4 

-1.65627950E-01 

3.78511216E-02 

3 

5 

-8.92789699E-02 

2,462(Ti301E-02 

3 

6 

-1.52776211E-02 

1.18819624E-02 

3 

7 

-4*77809758E-03 

4.74043228E-03 

3 

6 

-2.04304184E-03 

3.73013038E-03 

3 

9 

-1.30731458E-02 

9.95308345E-03 

3 

10 

-3.08981740E-03 

4.42260538E-03 

3 

11 

-9.07982028E-02 

8,r5^845866E-62 

3 

12 

3.20195257E-03 

5.04209210E-03 

3 

13 

•1.84017549E-03 

1.4129Q677E-02 

3 

14 

4.58966359E-03 

9.57498621E-03 

3 

15 

0. 

0« 

4 

1 

-3.88413133E-02 

-1.53173113E-02 


Appendix A 


Cont'd for Example 2 

PART 3-OERIVE CMATRIX 


LISTING OF POTENTIAL DISCONTINUITIES 


ROW 

COL 

REAL 

IMAGINARY 

1 

1 

7.19240047E-01 

-1.96046962E-01 

1 

2 

1.12190693E+00 

6. 19647631E-02 

1 

3 

7.82588972E-01 

4*45566363E>02 

1 

4 

l,11157379E+00 

2.77268333E-01 

1 

5 

7.49497603E-01 

1.90192592E-01 

1 

6 

9,29313913E-01 

3.40077119E-01 

1 

7 

5.99494642E-01 

2.43762016E-01 

1 

8 

6.58601669E-01 

3. 10765401E-01 

1 

9 

1.17338193E+Q0 

3* 57908726E-01 

1 

IQ 

9#74568491E-01 

4.05036350E-01 

1 

11 

1,32432438E+D0 

3.45992226E-01 

1 

12 

6.74430936E-01 

3.17596219E-01 

1 

13 

1.2l77d906E+0Q 

2.89990190E-01 

1 

14 

1.00525883E+00 

3.73063626E-01 

Z 

1 

2.28537254E-02 

-1.00678125E-02 

z 

2 

6.94409307E-02 

1.85305132E-02 

z 

3 

4.68362780E-02 

1.23600268E-02 

z 

4 

8.94395525E-02 

9.72468758E-02 

z 

5 

5.83969759E-02 

5.76615346E-02 

z 

6 

8.08251905E-02 

2* 00663205E-01 

z 

7 

5.21469123E-02 

1.11077066E-01 

z 

6 

5.64875102E-02 

1.93862280E-01 

2 

9 

9.92489587E-02 

4.33071999E-01 

2 

10 

7.97854832E-02 

4.95524806E-01 

2 

11 

1.12070872E-01 

1.21716905E-01 

2 

12 

6*84566594E-02 

2*20200364E-01 

2 

13 

1.26839334E-01 

5.48240810E-01 

2 

14 

1.15338995E-01 

5.66342821E-01 

3 

1 

-3.31905867E-03 

-2.68749724E-02 

3 

2 

4.05511239E-02 

-6.60891619E-02 

3 

3 

2.71556190E-02 

-4.47067358E-02 

3 

4 

1.26424457E-01 

-6.43118486E-02 

3 

5 

7.66394524E-02 

-4, 34688242E-02 

3 

6 

'2,27806'48iE-Oi 

-2,61098239E-02 

3 

7 

1.28759268E-01 

-2. 07230550E-02 

3 

8 

2.13937738E-01 

2, 25854757E-03 

3 

9 

4«66146099E»01 

-7,67512014E-03 

3 

10 

5.23098757E-01 

3.74847134E-02 

3 

11 

1.59097132E-01 

-7.80808271E-02 

3 

12 

2.44860925E-di 

7.19108775E-03 

>3 

13 

5.92013163E-01 

1. 00081713E-01 

3 

14 

6.07729750E-01 

8.19817326E-02 


Appendix A 
Cont’d for Example 2 


CHATRIX FOR N- 1 

NAIHH Ndr» •7dOHQ fRAs i.ood&d 


ROM COL real 

1 1 1.93944G33E’f00 

1 2 -3.S4698366E-02 

1 3 7.67700013E-01 

2 i 2.2T4l554&E-tii 

2 2 -1.27653711E-01 

2 3 9.33188544E-01 

3 1 8.43453971E-01 

3 2 -6.05133618E-02 

3 3 4.95997098E-01 

^ ^ 2 ^ ^ 

4 2 -1.89865173E-02 

4 3 2*73564386£-02 


IMAGINARY 

2.42293213E'f00 

7.71292667E-01 

3.13311215E-01 

’"i.32273833E-»-d0' 

5.68402451E*0I 

3.47877202E-01 

6.39851554E-01 

5.13395599E-01 

2.54582146E-01 

“5795289405E-02 

3.36392790£*02 

4.14206999E-02 
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Appendix A 

Lijstlng of Progreon PLOTT PI 


PRQGRAM FLUTT ( I NPUT i OUTPUT « TAPES=INPUTfTAPEfe=OUTPUT > 

DIMENSION DATA(600 ) « RH ( 10 ) « FRA ( 10 ) « NVECT ( 5 ) « VECT ( 5 ) 

REAL NKL 
COMPLEX A (6130) 

EQUIVALENCE (NDATAf DATA) 

COMMON/BLOCKA/NDATA(dOO ) t Af NRPf RMvFRAf DERSf OCt ALSQ«SQAL 

COMMON/BLOCKB/XKpLtSKOL.LSYritNpLANEtNfNKL 

COMMON/BLOCKC/NVECT rVECT 
CALL FTNBIN(0f 0«0) 

C * PROGRAM FOR UNSTEADY AIR FORCES USING * 

C ♦ OOWNWASH-VELOCITY POTENTIAL METHOD ♦ 

,C ♦..J.Y..J,K.v. HAVJJ-A of VA ♦. . . 

2 DO 3 J=1*6130 

3 A(J)S(0«<0.) 

DO 4 J=l»61 

4 NDATA(J)=0 

READ_L5i7 00) 

IF (EOF, 5) 999,5 

5 REAU(5»600) IC,NPDiNPWfNPA»NPP 
IF( IC.NE.l ) GO TO 999 
WRITE(6«fi2Q ) 

WRITE (6,700) 

WR,ITE(6.,fl00J. XCj.NPp t.NPWvNPArNPP 

REAiO( 5«80 0) IC, (NDATAl J) f J=2*56«4) 

WKITE( 6*800) IC* (NDATA( J) *Js2«56*4) 

IF ( IC,NE,2) GO TO 999 

REAO(5*eOO) IC, (NDATA( J) * J=3,59,4) 

WRITE (6,600) IC* (NDATA( J) * J=3«59*4) 

„IF( IC...NE,.3) _GP.X0 .999 

READ (5* 800 ) IC* (NDATA( J) * J=4*60 ,4 ) 

WHITE (6*800 I IC, ( NDAT A ( J ) * U=4 « 60 , 4 ) 

IF(IC,NE.4) GO TO 999 
READ(5f600) IC , NDMAT * NWMAT « NAMAT 
WRITE (6, 600) IC, NDMAT, NWMAT, NAMAT 

_ IF(IC*NE,5)_G0.. TO_ 999 

RE AD ( 5 * 60 0 ) IC , NMODEl t NM0DE2 * NKP * MTOT * NTE 
WRITE (6, 800) IC, NMODEl *NM0DE2«NRP* MTOT, NTE 
IF(IC.NE,6) GO TO 999 
NDATA(1)=62 
DO 25 J=l*15 

25. NDAT A ( 4* J+1 ) ^( NDATA ( 4»J-2 ) +NDAT A ( 4* J-1 ) ) ♦NDAT A ( 4* J ) + 

iNbATA('4*J-3) 

IF(NDAtA(61) .LE.800 ) 60 TO 30 

WRIJEJL6*105) 
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Appendix A 

Cont'd for Program FLUTT 

105 format (EX*LINE COUNT EXCEEDS STORAGE FOR DATA MATRIX*) 

GO TO 999 

30 IF (NPD .EQ.O) GO TO 35 

WRITE <6»106) (NDATA(K) *K=lf61) 

106 FORMATdX* DATA MATRIX ELEM£NTS*/(4X«4I4t 4m*4X*4I4) ) 
35 GO TO UiOf 120,130,I40),NDMAT 

110 CALL DMATR(MTOTiNPDfNTE) 

60 TO 500 

120 CALL DREAD (NPD) 

60 TO 500 
130 CALL DBLKl 

GO TO 951 
140 CALL DBLK2 
GO lU 951 

500 NM0DP-NM0DE1+NM0DE2 
NTOTW=NRP*NMOOE 
N1=NRP+1 
NADUW=N1*NKP 
NS12El=NKP*Nr'iODEl 
NADDW?=NADDW+MSIZEl 
IF(NMODE2.INJE,0) GO TO 510 
NTEMPrNMODEl 
NSIZE2=NSIZE1 
GO TO 520 
510 NTEMP=NM0DE2 

NSIZE2 = (\fRP+NM0DE2 
520 NAUDW3=NADDW2+NSIZE2 
NAD0WF = NADDW34-NSI2E2 
NCSIZE=NM0DE1*NTEMP 
NADDCrfJAUUWF + NSIZE2 
NT0T=rgAU0C+NCSI2E 
IF(NTOi .LT.llOOO) GO TO 10 
WRITE(6»e30) NTOT 
GO TO 999 

10 GO TO (210»220»230,240) vNWMAT 
210 CALL WHAT(NM00E«MT0T,NADDW«NT0TW) 

GO TO 310 

220 CALL WREAD(WMODE»NADDW*NTOTW) 

GO TO 310 
230 CALL WBLKl 
GO TO 953 
240 CALL WPLK2 
GO TO 953 

310 IF(l\»PW.tO»0> GO TO 45 
WRITE(6,840 ) 

LALL ALIST (AsNADDW.NMODEtNRP) 

45 READ(5*810) IC^NC. (RM«I) *I=ltNAMAT) 

WRITE(fa*810)ICtNC» (RM(I) fI=itNAMAT) 

IF(IC,NE.29) GO TO 999 

READ (5* 610) IC«NCf (FRA (I) ♦ 1=1 f NAM AT) 

WRITE(6t0lO) IC*NCi (FRA(I) tI=l,NAMAT) 
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Cont'd for Program FLUTT 


IFdC.NE.iO) GO rO 999 
DO 940 N=1*NAMAT 
WRITE (6»8S0) 

WRITE (6*860) N*NRP*RM(N) *FRA(N) 

DO bO J=1*NADDW 
C CALL ACALC 

5Ci A(J) = (0,*U.) 

CALL ACALC(MTOT) 

IF(NPA.NE.O) CALL ALJ ST ( A * 0 * NRP * Ml ) 

CALL cMATR(Nl*NT0TW*NAD0W»NP10DEl*NM0DE2fNTEMP*NSIZEl* 
lNSIZE2*NCSIZE*NADDC*NPP*NADDW2tNAD0W3*NADDWF) 

700 FORMAT UX*79H 
1 

800 FORMAT ( IX * 12*2013) 

810 FORMAT ( IX * 12* 4X * I3*6F10.5) 

820 format ( ?X//1X*PK06RAM FOR UNSTEADY AIR FORCES USING*/ 
12X* dOWNWASM-VELOOiTY POTENTIAL METHOD*// 

21X*BY J*K. HAVILAND AND Y.S, YOO • UNIV OF VA*// 

630 FORMAT (1X*ST0RAGE CAPACITY FOR TOTAL A MATRIX lS*tl6) 
640 FQKMATCIX//* WMATRIX. ELEMENTS*// ) 

850 FORMAT (2X//1X* PART 2 - DERIVE AERODYNAMIC MATRIX*/) 
860 F0RMAT(1X*A MATRIX**I3** NRP=**I3** MACH NO=**F8*5»* 
1FRA=*,F8.5) 

940 CONTINUE 

950 CONTINUE 

951 CONTINUE 
953 CONTINUE 

GO TO 2 
999 CONTINUE 
STOP 
END. 


SUBROUTINE DMATR ( MTOT « NPD * NTE ) 

DIMENSION UATA(800 ) * XO ( 4 ) * YO ( 4 > * ZO ( 4 ) * COSG ( 
XDELX(14) 

DIMENSION RM( 10) *FRA(10) 

COMPLEX A(6130) 

DIMENSION NSUR<4) *NSYMZ(4) *NSYMY(4) *NGE0M(4) 
DIMENSION SFORU4) *SAFT(14) , SIN (14) *S0UT(14) 
DIMENSION iMb(14) ♦ NOT ( 14 ) * NSURF ( 14 ) 
EQUIVALENCE (NDATA*DATA) 

COMMON/BLOCKA/NDATA(800 ) *A,NRP*KM*FRA*DERS* 

DO. 950 J=l* 15 

NADD=N0ATA(4*J-3) 

NFIX=MDATA(4*J-2) 

-NFL=NdATA(4*J-1) 


4) *SING(4) ( 


«DELSU4) 

0C« ALSQtSQAL 


NL=NDATA(4*J) 

IF(NL.EO.O) GO TO 950 

GOTO (110 *115 *120 *120 *120 *120 *12 0*120 *120 *120. 120*120* 

1120*120.120) ,J 
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X.XO RE:AP(5».6P0) IC* (NSUR(I)tI=l*NU 
WRITE (6t 800) ICf (NSUR(I) f 
IFUC.NE. 7) STOP 
READ(5»800) iLi (WSYMZ(I) .I=lfML) 
WRITE (6» 800 ) IC* ( NSYM7 ( I ) t I=X t NU 
IF{IC.NE, 8) STOP 
READ(5t600) IC * ( NSYMY ( I ) • 1=1 tNL) 
WRITE(fa«fiOO)ICi (NSYMY( I ) i 1=1 «NL) 
IFdC.NE, 9) STOP 
READ(5»e00) ICt (NGE0M(I)<I=1*NL) 
WRITE (6«600)IC* (NGE0M( I ) f I=lfNL> 
IF(IC.NE.IO) STOP 
READ 610) IC,NCt (X0( I ) * I=lf NU) 
WRITE(6t810)IC*NCi (XO(I) ♦I=l»NU 
IFdC.NE.ll) STOP 
READ(5t610) iCtNCf (YOU) ♦I = ltNL) 
WRITE(6t610)ICfNC« ( YO ( I ) « 1=1 « NL ) 
IF(IC,NE,12) STOP 
READ(^f610) ICtNCt (ZOU ) • I = 1«NU 
WRITE(6.fli0)ICiNC* (Z0<I) .I=liNL) 
IF(IC.NE.13) STOP 

READ< 5*610) IC,NC» (COSG(I) *I=1.NL) 
WRITE(6*610 ) ICfNC* (COSG( I ) * 1 = 1* ML ) 
IF(IC.NE,14) STOP 

KEAD(5*610) ICtNC* (SING( I) *I=1*NL) 
WR1TE(G*810)IC*NC* (SlNGd ) <I = 1*NL) 
IF(IC,NE,15) STOP 
GO TO 500 

115 READ(5*800) IC « ( NSURF ( I ) * I = 1 * MTOT ) 
WKITE(6»800 ) IC* (NSURF< I ) ♦ I = 1*MT0T ) 
IF(IC.NE*16) STOP 

REAU(5*800) IC* (NQT( I ) *I=l*MTOT) 

WRITE (6*600 )IC* (NOT (I ) ♦ I=l*MTOT) 

IF(IC,NE.17) STOP 

READ(5*800) IC * ( NB ( I ) • 1 = 1 * MTOT ) 

WRITE(6*800) IC* (NB( I ) * 1 = 1* MTOT) 

IF(1C.NE.18) STOP 

ND=0 

DO 10 I1=1*MT0T*6 
12=11+5 

IF(12*GT,MT0T) I2=MT0T 

READ (5* 81 0 ) ICf NC* (SFOR{ I ) f I=Ilf 12) 

WRITL(6ffllO) ICfNCf ( SFOR ( I ) , l=ll f 12 ) 

ND=NU+1 

10 IF( IC,NE,19,0R,NC.NE,ND) STOP 
ND=0 

DO 20 I1=1*MT0T*6 
12=11+5 

IF(I2.GT.MT0T) I2=MT0T 

RE AD( 5*810) ICfNC* ( SAFT (I ) * I = I 1 * 1 2 ) 

WRITE (6* 810) IC*NC, (SAFT(I)f 1=11*12) 
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NOsND<fl 

20 IF(IC«NC.20«OR.NC.NE*NOI STOP 

NOsO 

DO 30 llsi,nT0Tt6 

I2SI145 

IF(I2.6T.nT0TI X2snTOT 
RCAOISfSlO) ZCfNC«(SlN m«Z«XltX2) 

URITC(6«aiO) XC«NC«(SXN a)fIaZl«Z2) 

_ NDsND+1 

30 IFnC*NC.21*0R.NC*NC*ND> STOP 
NOsO 

00 40 Xlsl,nT0Tt6 
I2S1145 

XF(I2»6T.nT0T) X2SRT0T 

RCAOfStdlO) ICtNC* ISoUTIDf Xall«12| _ 

WRlTCfStSlO) IC«NCf (S0UT(X)«ZsIlfI2| 

N0sND4l 

40 IF(IC.NE.22.0R*NC*NE»N0) STOP 
0ERSsS0UTm»SZN(ll 
IF(XC»NE.22) STOP 
120 GO TO 500 
500 DO 940 NLINEsl*NL 

K IsN A00> i NF I X^NFL ) « ( NL INE*X I 

K3SK1+NFXX 

K2SK3-1 

K4SK2+NFL 

GO T0(210*220f 230*240 « 250*260 f 250 t 260 f 250 t 260*250 #240 f 

1330*340«350) * J ” 

210 NOATA(K1)sNSUR(NLZNEI 

NOATA ( K14-1 ) sNSYRZ ( NLZNE ) 

N0ATA(K1^2)3NSYNY INLINE) 

NOATA ( K2 ) sNGEOn ( NL XNE ) 

0ATA(K3)=X0(NLINE) 

DATA I K3+1)=Y0 INLINE) 

DATA(K3^2)=Z0INLINE) 

DATA I K34-3 ) rCOSG INLINE) 

DATA I K4)sSING INLINE) 

60 TO 32 

220 LINEsQ _ _ 

DO 930 «sl*MTOT ' “ ' ‘ “ 

NQsNQTIR) 

DO 920 ZQsi*NQ 

IFII0.LT«NQ«0R .NBIPU.LT*!) go to 222 
lF(NLINt.LE.INRP«NTE)) GO TO 930 
LINEsNLINE 

ntb*o 

do 4 MTESTsl*MTOT 
ZFlNB<nrEST).GE.l) NTBsNTB<fl 
XFINTB.EQ. INLINE->NRP4NTE) )nsNTEST 
IF{NTB*EQ.INLINE«NRP4NTE)>G0 TO 224 
4 CONTINUE 
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222 LINE=LlWt: + l 

IF(LlNE*r\.E,NLINE) GO TO 920 
224 DELX(M)=(SAFTJM)-SF0R(M) )/NQT{M) 
DELS(n)=SOUT(~M)-SIN(M) 

ITYPE=J-1 

60 T0(42 0 f 430«440<450<460 t4SOt460«450f460« 
1450 » 460 » 530 » 540 *550 ) ♦ I TYPE 


230 

GO 

TO 

220 


240 

GO 

TO 

260 


250 

GO 

TO 

220 


260 

NLINE; 

= NRP- 

-NTE+NLINE 


GO 

TO 

220 


330 

GO 

TO 

950 


340 

GO 

TO 

950 


350. 

GO 

TO 

220 


420 

NDATA 

(Kl): 

:LINE 


NDATA. 

(K2) = 

:NSURF.(F1) 


DATA(K3)=SFOR(M)+( ,75+ (lO-l) )*DELX(W) 
DATA (K4 )=SIN(M)+PEUS(M)*,5 
GO TO 32 
430 GO JO 450 
4*tu GO TO 4bU 
450 NljATA(Kl)=LINE 

XF=SF0R(M)+0ELX(M)+,25 
DATA(K3)=XF+( IQ-X)^DELX(M) 
DATA(k3+1)=OATA(K3)+DELX(M) 

DATA(K3 + 2)=SIN|W). 

DATA(K4)=S0UT(M) 


IF(IQ,_EQ,iNJQ,AND.NB.(M),LE,1) go to 452 
GO TO 453 

452 DO 6 nT=l»MTOT 

lF{SlN(nT) ,LT,DATA(K4) . AMD • SOUT ( FIT ), GT , DATA ( K3+2 ) 
l,AND.riT.GT,M) GO TO 5 
GO TO 6 

5 DATA(K3+1 )=$FOR(MT)+(SAFT(MT)-SFOR(MT) )/(N0T(MT)*4) 
GO TO 453 

6 CONTINUE 

453 GO TO 32 

460 NDATAfKl )=LINE 
NDATA(K2)=NB(M) 

DATA ( K3 ) =SAFT (M) -DEUX (M)*. 75 
DATA(K3+1 )=DATA(K3)+DELX(M)*NDATA(K2) 

DATA(K3+2 )=SIM(M) 

DATA(K4)=S0UT(M) 

NL1NE=NLINE-NRP+NTE 
GO TO 32 
530 GO TO 950 
540 GO TO 950 
550 GO TO 450 

32 IF(iMPD,NE.O) WRITE(6,107) J « (NOATA(K) *K=KlfK2l 
107 F0RMAT(3Xt7l4) 
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IF(NPD,NE.0) WRITE(6tX09) J t ( DAT A ( K ) f K=K3 , K*f ) 
109 F0RHAT(2Xf I3f 5F13.8) 

600 FORMAT(lXtI2»20I3) 

610 FORMAT ( 1 X « 12 « 4X 1 13 1 6 F 10 . 

GO TO 9"46 
920 COWTlTaUE 
930 CONTINUE 
940 CONTINUE 
950 CONTINUE 

RETURN 

END 


subroutine DKEAD(NPD) 

DIMENSION DATA(600 ) « RM ( 10 ) f FRA ( 10 ) 

COMPLEX A (61301 
EQUIVALENCE (NDATA* DATA) 

COMMON/eLOCKA/NDATA(600 ) t A t NRP ♦ RM . FRA t UERS . OC t ALSO* SOAL 

35 DO 5 J=62*800 
5 NDATA(J)=0 

ND1=6 

ND2=7 

DO 50 J=l*15 
NAUD=MDATA(4*J-3) 

NFIX=NDATA(4*J-2) 

NFL=N0ATA(4*J-1) 

NL=NDATA (4*J) 

IF(NL.EO.O) GO TO. 50 

ND1=NU1+1 

ND2=ND2-H 

DO 40 NLINE=liNL 

K1=NADD+(NFIX+NFL)*(NLINE-1) 

K3=Kl+NFIX 

K2=K3-1 

K4=K2+NFL 

READ(5*800) IC i I T t ( NDATA ( K ) t K=K1 ♦ K2 ) 

IF( IC.EQ. (NDl } .AND.IT^EQ.J) GO TO 37 
WRITE(6«106) J«IT 

106 FOR«AT( lX*TYPE*t I3**CARD OUT OF ORDER READ FIXE0«tI3) 

GO TO 999 

37 IF(NPD.NE.O) WRITE(6«d00) IC f I T t ( NDATA ( K ) • K=K1 t K2 ) 

READ (5« 810) ICf IT« ( DATA ( K ) • K=K3 « K4 ) 
DERS=DATA(K4)-DATA(K3+2) 

IF(IC.EQ.(NU2 ) •AND. IT.EQ* J) GO TO 36 
WRITE(6«110) JfIT 

110 F0RMAT(lX*TYPE*tI3«*CAR0 OUT OF ORDER READ FL0AT*tI3) 

60 TO 999 

36 IF(NPD.NE.O) WRTTE(6.610) IC , IT » ( DATA ( K ) • K=K3t K4 ) 

40 CONTINUE 

ND1=ND1+1 

N02=ND2+1 
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30 CONTINUE 

800 FOKMAT (lX*I2t20I3) 

6X0 FORMAT(lX*I2«4XfI3»6F10.5) 
999 CONTINUE 

return 

END 


SUBROgj If'JE WMAT ( NMOOE i NTOT t NAODW » NTOTW ) 

DIMENSION OATA(800 )f NSUR (14 ) « SIN U4 ) f SOUT ( I<f ) « XROT< 8«f ) 
DIMENSION RM(10» fFRA(lO) »H(84) ,ALPH(84) 

COMPLEX A(6130) 

EQUIVALENCE (NDATA» DATA) 

COMMON/BLOCKA/NOATA(800 ) • A«NRPtRMt FRA*OERS«OCt ALSQtSQAL 
DO 5 J=1»NT0TW 
5 A(NAOCW-kJ) = (0) 

NMOMT=NMODE*MTOT 

N0=0 

DO 23 lX=ltNM0MTf6 
12=11+5 

IF(I2,GT,NM0MT) I2=NM0MT 
READ(5«B10) IC«NC«(H(l)fI=Il<I2) 

WRITE (8« 810) ICfNCf (H (I) t 1=11 i 12) 

NO=ND+l 

23 1F(IC,NE.23.0R.NC.NE,NO) STOP 
NU=0 

DO 24 Zl=lfNM0MTf6 
12=11+5 

1F(I2.GT.NM0MT) I2=NM0MT 
RLAD(5f810) IC,NC* (ALPH(I),I=I1,I2) 

WRiTE(6«8lO) ICtNCf (ALPH(I) f I=I1« 12) 

N0=ND+1 

24 IF(IC,.NE.24.0R.NC.NE.ND) STOP 
Nu=u 

DO 25 Il=lfNM0MT»6 
12=11+5 

IF(I2.GT,NM0MT)I2=NM0MT 

READ (5« 610) ICfNCf (XROT( I) . I=tlf 12) 

WRITE (6«810) IC«NC« (XRpT(I) «IsIl«I2) 

ND=NU+1 

25 IF(IC.NE.25.0R.NC.NE.ND.) STOP 
ND=0 

00 26 Il=ltMT0T»6 
12=11+5 

1F( 12.GT.MT0T) 12=MT0T 

REAU(5f 610 ) ICiNCt (SIN( I ) « I=Il« 12) 

WR1TE(6«610) IC«NCi(SlN (I),IsIl«I2) 

ND=ND+1 

26 IF(1C,NE.26.0R.NC.NE.ND) STOP 
ND=0 
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D0_ 27 = 

I2=li+5 

IF< I2.GT,MT0T) I2=MT0T 
•^LAU{St810) IC,NC» (SOUT(I) f I = IltI2) 

WRlTE(6i610) IC«NC«(S0UT(I)«I=11«I2) 

NU=ND+1 

27 IF( IC,NE,27.0R.NC.NE,ND) STOP 

KEAO(5»800) IC# (NSURt 1 ) f IslfMTOT) 

WRITE(6»800) IC* (NSUR( I ) * I-l»MTOT) 

IFaC,NE.28) STOP 
800 FORMATdXf I2«20I3) 

610 FOHriAT( IX f I3«6F10,5) 

pO__ qOJ'IODErl .NPIO^ 

DO 40 LIN2=lfNRP 

DO 40 M=lfMTOT 

NADD2 = NDATA(5)+<f*(LIN2-l) 

IFANSURtMl. *Ni:,.J4QMAlNADD2+l).)._^ TO 40 

X= DATA(NADD2+2) 

S1=DATA (NADD2+3) 

IF(SlPJ(M) .GT.Sl .OR,(SOUT(M) .LT.SD) GO TO 40 

N=MTOT*(MODE-l)+M 

J=NRP* ( mode-1 ) +L IN2+NA0DW 

A ( J ) =CMPLX (_< ALPH ( NJ ). t ( H ( N ) + ( X-XROT ( N ) ) »ALPH ( N ) ) ) 

40 CONTINUE 
999 . CONTINUE 
RETURN 
ENU 


SUBROUTINE WREAD ( NMODE f NAOOW t NTOTW ) 

COMPLEX A<G130) 

DIMENSION OATA(800 ) »RM ( 10 ) t FRA ( 10 ) 

EQUIVALENCE (NOATA I DATA) 

COMMON/BLOCKA/NOATACeOO ) *A«NRP«RMt FRA«DERSf OCt ALSQtSQAL 
DO 5 J=l#NTOTW 
5 A(NADDW4J)=(0> 

JsNRP*NMOOE-i>NAODW 

ND=0 

IWsNADOW>d 
DO 10 Il7lWfJ«3 
I2SI1+2 

IF(I2.GT.J) I2sJ 

READ(5<810) ICfNCf (A(I) f IsIl«I2) 

WRITE (6« 810 )ICfNC« ( A ( 1 > t I^ll* 12 ) 

ND=ND4l 

IF(IC.NE.23.0R.NC,NE.ND) STOP 
10 CONTINUE 

810 FORMATUX*I2*4Xf I3f6F10,5) 

RETURN 

END 
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SUBROUTINE ACALC 

DIMENSION DATA(800 ) • RM ( 10 ) f FRA ( 10 ) « NVECT ( 5 ) • VECT ( 5) 

REAL NKL«NLK«NKLNLK 
COMPLEX A(6130) f AK,AA 
EQUI VALENCE (NDATA* DATA) 

COMMON/BLOCKA/NDATA(800 ) iA,NRP»RMf FRAfDERSf OCf ALSOfSOAL 
COMMON/BLOCKB/XKOL«SKOLtLSYM»NPLANEtNfNKL 
COMMON/BLOCKC/NVECTfVECT 
ALSQ=1.0-RM(N)*RM(N) 

SQALsSQRT(ABS(ALSQ) ) 

OC=RM(N)*FRA(N) 

NL1=NDATAU) 

NL2=N0ATA(8) 

DO 900 LlNlK=lfNLl 
NADD1K=NDATA(1)+9»(LIN1K-1) 

NSURK=NDATA(NADD1K) 

XOK=DATA(NAOD1K44) 

Y0K=DATA(NADDlK+5) 

Z0K=DATA(NADD1K+6) 

COSGK=DATA(NADDlK+7) 

SINGK=0ATA(NADD1K+8) 

DO 900 LIN2K=XfNL2 
NAQD2=N0ATA<5)+‘f*(LIN2K-l) 

IF (NOATA (NAD02+1) tNE.NSURK) GO TO 900 
KsN0ATA(NAD02) 

XK=0ATA{NADD2+2) 

SK=0ATA(NADD24-3) 

00 898 LlNlLsl , NLl 
NA0D1L=NDATA(1)+9*(LIN1L-1) 

NSURL=N0ATA(NADD1L) 

NSYMZ=NDATA(NADD1L+1) 

NSYMYrNDATA(NADDlL-f2) 

NGE0M=NDATA(NA0D1L+3) 

XOL=DATA(NAODlL+4) 

Y0L=DATA(NADD1L+5) 

Z0L=DATA(NADDlL+6) 

C0SGL=DAT A ( NAODlL+7 ) 

SINGL=DATA(NA0DlL+6) 

DO 897 MSYM=1»4 

GO TO (105f 110tll5«120) fMSYM 

105 KSYMZsKSYMYaLSYMZsLSYMYrl 
GO TO 130 

110 IF(NSYM2,NE,0) GO TO 111 
MSYMsS 

60 TO 897 

111 KSYMZa-l 
LSYMZsNSYMZ 
GO TO 130 

115 IF(NSYMY.EQtO) GO TO 896 

116 KSYMYs-1 
LSYMYsNSYMY 
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GO TO 130 
lao KSYHZsLSYMZai 
GO TO 116 

130 LSYM=LSYMY*LSYP1Z*KSYMZ 
NPLANEsKSYMZ*KSYMY 
IF ( NSURL • NE . NSURK ) NPLANE=0 

IF ( ( ABS ( SI NGL ) ♦ < 1 -KS YMZ ) +ABS ( COSGL ) * ( 1-KSyMY ) ) . NE • 0 ) 
inplane=o 
XK0L=XK+X0K-X0L 
delz=zok»ksyhy-zol 

OELY=YOK*KSYMZ-YoL 

T1KLsC0SGK*C0SGL»KSYMZ+SINGK*SINGL*KSYMY 

skol=delz*singl+dely+cosgl+sk*tikl 

SINGKLsSINGK»COSGL*KSYMY-COSGK*SINGL+KSYMZ 
NKL=DELZ*C0SGL^DELY*SINGL+SK^SINGKL 
DO 897 LIN2L=1»NL2 
NADD2=NDATA(5) + »t*(LIN2L-l) 

IF(NDATA(NADD2+1) .NE. NSURL) GO TO 897 
LsNDATA(NADD2) 

XL=DATA(NA0D2-*-2) 

SLsDATA(NADD2>3) 

IF (NPLANE.NE.O ) GO TO 153 

NLK2-SL^SINGkL-0ELZ»C0SGK*KSYMZ+DELY*SINGk*KSYHY 
NKLNLKsNKL+NLK . 

153 DO 897 NTYPEs3il5 

call TEST(NTYPE.L.NL*NADDJ) 

IF(NL.EO.O) GO TO 897 
IF<N0ATA(NA0DJ) .NE.L) go to 897 
JTYPE=NTYPE-2 

GO T0(230f 240«250.260«270f 280. 290 .300* 3lO« 320 .330.340 
1.350) • JTYPE 

230 call AC0N3(T1KL.AA .NKLNLK .AK.XL,SL.L.K,NADDJ) 


GO TO 896 

240 CALL AC0N4(T1KL.AA . NKLNLK . AK . XL • SL. L . K . NADDJ ) 

GO TO 896 

250 CALL AC0N5(TiKL.AA . NKLNLK. AK.XL.SL.L.K.NaDOJ) 

GO TO 896 

260 CALL AC0N6(T1KL«AA .NKLNLK. AK.XL.SL.L.K.NaDOJ) 

GO TO 896 

270 CALL AC0N7(T1KL.AA .NKLNLK . AK« XL.SL.L.K.NAOOJ) 

GO TO 896 

280 call ACON8(TiKL.AA . NKLNLK. AK.XL.SL.L.K.NaDDJ) 

GO TO 896 

290 call ABLK9 
GO TO 897 

300 CALL ABLKIO 
GO TO 897 

310 CALL ACONll (TlKL.AA .NKLNLK. AK.XL.SL.L.K.NADDJ) 
GO TO 896 

320 CALL AC0N12 (TlKL.AA .NKLNLK. AK.XL.SL.L.K.NADDJ) 
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/ 


330 

340 

350 

896 


897 

698 

900 


GO TO 896 
CALL ABLK13 


GO TO 897 
CALL ABLK14 
GO TO 897 
call AC0N15 
CONTINUE 


CTlKLf AA fNKLNLKtAKfXLtSLtLtKrNAOOJ) 


NADO=(L+(K«*l)*(NRP-H) ) 

A(NAD0)sA(NA00)4>AA 

CONTINUE 

CONTINUE 

CONTINUE 


RETURN 

END 


SUBROUTINE AC0N5 (TlKLtAA tNKLNLK , AK*XL»SLfL»Kf NADDJ) 
DIMENSION DATA(800 ) t RM ( 10 ) • FRA ( 10 ) t VECT ( 5 ) t NVECT ( 5 ) 

REAL NKL*NlK*NKLNLK 
COMPLEX A(6130) f AKiAA 
EQUI VALENCE <NDATA» DATA) 

COMMON/BLOCKA/NDATA(600 ) t A t NRP * RM » FR A t DERS f OC f ALSO « SQAL 
COMMON/BLOCKB/XKOL,SKOLtLSYMiNPLANE»N,NKL 
COMMON/BLOCKC/NVECTfVECT 
RU*S) = (XK0L-X)**2^ALSQ*( (SK0L-S)**2+NKL**2) 

AA=0 

IF(NDATA(NADDJ) tEOtL) GO TO 100 
GO TO 894 
100 XFOR=VECT(l) 

XFORW=XFOR 

XLSTAR=XL 

XAFT=VECT(2) 

SIN=VECT(3) 

S0UT=VECT(4) 

ALPs( (XAFT-XFOR)*(SOUT-SIN) )/12,566371 

DELXsO 

NBOXES=l 

GO TO 267 

ENTRY AC0N6 

AA=0 

IF(NDATA<NAOOJ) .EQ.L) GO TO 200 
GO TO 894 

200 NB0XESsNVECT(2) 

XF0KWsVECT(2) 

XAFTW=VECT(3) 

DELX=(XAFTW-XFORW)/NBOXES 

XFORsXFOKW-DELX 

XLSTAR=XF0RW*.5*0ELX 

XAFT=XFORW 

SINsVECT(4) 
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S0UT=VECT15) 

ALPsDCLX«<S0UT-SIN)/12, 566371 
267 DO 694 NBOX=l< NBOXES 
XFORsXFOR+DELX 
XLSTAHsXLSTAR-fOELX 
XAFT=XAFT+DELX 
RKLSQ=R(XLSTAR«SL) 

RKLsSQRT(RKLSQ) 

RKLCU=RKLSQ*RKL 
IF(NPLANE.NE.O) GO TO 440 
IF(L*EQ*K) 60 TO 450 
AKs*ALSQ«ALP/RKLCU 
GO TO 450 

440 AKsFUNK ( X aft « SOUT ) -FUNK < XAFT t SIN ) -FUNK ( XFOR « SOUT ) 4>FUNK 
KXFORfSIN) 

AKsaK«(1*/12» 566371) 

450 IF(OC«NE,0) AKsAK* < 14> ( 0. 1 1 • ) «OC«tRKL/AL$Q ) 

AKrAK*TlKL 
IF(L.EQ,K) GO TO 500 

IFCNPLANE.EQ.O) AKsAK<i-NKLNLK*ALP* ( -3«*ALSQ«*2/RKLSQ 
l-3.«(0. <!• >*0C«ALSQ/RKL+0C*«2)/RKLCU 
500 AKsAK^LSYM 

IF (OC.EO.O) GO TO 495 

XFUN=FRA(N)*(XL-XLSTAR-MRM(N)**2/ALSQM'(XK0L*XLSTAR)- 
1 (RM(N)/ALSQ)*RKL) 

AK = AK»CEXP( (0* »1. )*XFUN) 

495 AA=AA+AK 
694 CONTINUE 
RETURN 
ENU 


subroutine AC0N7 (T1KL«AA «NKLNLK« AKi XUf SL<L«Kf NAOOj) 
DIMENSION DATA<6Q0 ) « RM ( 10 ) • FRA ( 1 0 ) f VECT < 5 > f NVECT ( 5 ) 

REAL NKLtNLKiNKLNLK 
COMPLEX A(6130) «AK«AA 
EQUIVALENCE ( NDATA t DATA ) 

COMMON/BLOCKA/NDATA(600 ) i At NRP«RM«FRAf DERStOC* ALSQf SOAL 
COMMON/BLOCKB/XKoLtSKOLtLSYM^NPLANErtNtNKL 
COMMON/BLOCKC/NVECTtVECT 
R ( X « S ) s ( XKOL-X } *«2’*-ALSQ« ( < SKOL-S ) *4c2^NKL>»«2 ) 

AAsO 

IF(NDATA<NADDJ) .EQ.L) GO TO 100 
GO TO 894 
100 XF0R=VECT(1) 

XFORW=XFOR 

XLSTARsXL 

XAFT=VECT(2) 

SIN=VECT(3) 

S0UTsVECT(4) 
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ALPsMXAFT-XFOR)*(SOUT-SIN) )/12.566371 

OELXsO 

NUOXCSsl 

GO TO 267 

ENTRY AC0N8 

AASQ 

IF(NDATA(NADDJ) .EQtL) -GO TO 200 
GO TO 894 

200 NB0XES-NVECT(2) 

XF0RW=VECT(2) 

XAFTW=VECT(3) 

DELX=(XAFTW-XF0RW)/NB0XES 

XFOR=XFOftW-DELX 

XLSTAR=XF0RW-*5*DELX 

XAFT=XFORW 

SINsVECT(4) 

S0UT=VECT(5) 

ALP=DELX*(S0UT-SIN)/X2. 566371 
267 UO 694 NBOXsi « NBOXES 
,XFOR=XFOR-»-DELX 
XLSTAR=XLSTAR+DELX 
XAFT^XAFT+DELX 
RKLSQ=R(XLSTARfSL) 

RKLsSQRT(RKLSQ) 

RKLCU=RKLSQ*RKL 
IF(NPLANE,NE.O) GO TO 420 

ak=-also*alp/rklcu 

AKUs-OC*ALP/RKLSQ 

AKs < AK+ ( 0 • f 1 . ) *AKU ) ♦TlKt+NKLNLK*ALP* ( *3. *ALSQ**2/RKLSQ 
1 -3.+(0. »1. >*OC»ALSQ/RKL+OC**2)/RKLCU 
420 IF(NPLANE,EQ.O) GO TO 425 

IF<NPLANE.NE.O> AK=FUNK ( XAFT • SOUT ) -FUNK ( XAFTf SIN) •FUNK 
l(XFOR*SOUT)+FUNK(XFORfSlN) 

AK=AK*(1,/12. 566371) 

AKU=FUNKU(XAFT»SOUT)-FUNKU(XAFT«SlN)-FUNKU(XFOR»SOUT) 

l+FUNKU(XFORfSIN) 

AKU=AKU*(0C«RM(N)/<12«566371«SQAL) ) 
AKU=AKU+(OC/ALSQ)*(RW(N)*(XLSTAR-XKOL)+RKL)*AK 
AKs(AK+(0. )*AKU)*T1KL 

425 XFUN=FRA(N)« (XL-XLSTAR’('(Rn(N)««2/ALSO)«(XKOL•XLSTAR)• 

l (RM(N)/ALSQ)>ifRKU) 

AK=AK*LSYM+CEXPnO. *1. )♦ XFUN) 

AA*AA-t-AK 
694 CONTINUE 
RETURN 
END 
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SUBROUTINE ACONlKTlKLt AA tNKLNLKt AKtXLtSL«LiK«NADOJ) 
DIMENSION DATA(600 ) • RM ( 10 ) i FRA ( 10 ) t NVECT ( 5) « VECT ( 5) 

REAL NKLfNLKfNKLNLK 
COMPLEX AC6130) «AK«AA 
EQUI VALENCE (NDATAf DATA) 

COMMON/BLOCKA/NOATA(BOO ) «A«NRP«RM*FRAf OERSf OCtALSQfSQAL 
COMMON/BLOCKB/XKoLfSKOLtLSYMfNPLANEfNiNKL 
COMMON/BLOCKC/NVECT»VECT 
R ( X « S ) = ( XK0L*X ) ♦♦2'l>ALSQ4‘ ( ( SKOL-S ) «*2-l'NKL«*2 ) 

AA=0 

IF(NOATA(NADDJ)»EQ.L) 60 TO 100 
GO TO 894 
100 XF0R=VECT<1) 

XFOKWsXFOR 

XLSTAR=XL 

XAFT=VECT<2) 

SIN=VECT(3) 

S0UT=VECT(4) 

ALP=( <XAFT-XF0R)*(S0UT-SIN) )/12,566371 

DELX=0 

NBOXESsl 

GO TO 267 

ENTRY AC0N12 

AA=0 

IF(NOATA(NADDJ).Cfl.L) GO TO 200 
GO TO 894 

200 N80XES=NVECT(2) 

XF0RW=VECT(2) 

XAFTWsVECT(3) 

DELX=( XAFTW-XFORW)/NBOXES 

XFUR=XFORW-DELX 

XLSTAR=XF0RW-#5*DELX 

XAFTsXFORW 

SINsVECT(4) 

S0UT=VECT<5) 

ALP=DELX*<SOUT-SIN)/12. 566371 
267 DO 894 NBOX=l f NBOXES 
XFORsXFOR+DELX 
XLSTAR=XLSTAR+DELX 
XAFT=XAFT+OELX 
RKLSQ=R(XLSTAR«SL) 

RKLsSQRT(RKLSQ) 

RKLCUsRKLSQ*RKL 

AK= U * /12 • 56637 ) * ( FUNK ( X AFT t SOUT ) «FUNK f XAFT t SIN ) -FUNK 
1 ( XFOR « SOUT ) 4>FUNK < XFOR t SIN n 
AKUs ( FUNKU ( X aft f SOUT ) -FUNKU ( XAFT # SI N ) -FUNKU ( XFOR f SOUT ) 
I'i'FUNKUtXFORtSIN) )«(FRA(N)/(12.566371*SQAL) ) 

AKU=AKU+ (RKL*RM(N)+XLSTAR-XKOL)*AK*<FRA(N)/ALSO> 

AKs(AK+(0. *1. MAKUli^TlKL 
IFCNPLANE.I^.O) go to 425 
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1F(L«E;0.K) go to H25 

AKs AK +NKLNLK*ALP*(-3.*ALS0**2/RKLS0 

1 )*0C*AtSQ/RKL+0C**2)/RKLCU 

^25 XFUN=FRA(N)*(XL-XLSTAR+(RM<N)**2/ALS0> 
l*(XK0L-XLSTAR)-(RM«N)/ALSO>*RKL) 

AK=AK*LSYM*CEXP(<0.»1#)* XFUN» 

AAsAA-i-AK 
694 CONTINUE 
RETURN 
END 


FUNCTION FUNKCXfS) 

REAL NKLtNLKfNKLNLK 

DIMENSION OATA(600 ) tRMdO ) «FRA< 10) «NVECT(5) t VECT(5) 
COMPLEX A(6130) 

EQUIVALENCE (NDATAf DATA) 

C0MMON/BLOCKA/NDATA(600 ) t AfNRPt RMt FRAi DERS«OC« ALSQtSQAL 
COMMON/BLOCKB/XKOLfSKOLfLSYMtNPLANE»N*NKL 
R«XtS)=(XK0L-X)»*2+ALSQ*( <SK0L*S) ♦♦a+NKL^^a) 
IF(NPLANE.NE*0) GO TO 10 
IF(ABS(NKL/DERS)*6Tf0001) 60 TO 5 
GO TO 10 

5 FUNKS (1,/NKL)«ATAN2(NKL*SQRT(R(X«S) ) « < (XK0L-X)«(SK0L-S) ) ) 
RETURN 
10 CONTINUE 

IF<ABSUXKOL-X)/DERS)«LT«*0001) GO TO 20 
IF(ABS( (SK0L-S)/DERS) •LT«*0001) 60 TO 20 
FUNKsSQRT(R(XfS) )/( (XK0L»X)«(SK0L->S) ) 

60 TO 30 
20 FUNKsO 
30 CONTINUE 
RETURN 
END 


FUNCTION FUNKU(X*S> 

DIMENSION DATA(800 ) «RM( 10 ) «FRA( 10 ) t NVECT(5) *VECT(5) 
COMPLEX A(6130) 

EQUIVALENCE <NOATA« DATA) 

REAL NKLfNLKtNKLNLK 

COMMON/BLOCKA/NDATA(600 ) «A«NRPfRM«FRAfDERSfOC*ALSQfSQAL 
COMMON/BLOCKB/XKOLtSKOL«LSYMfNPLANE*NfNKL 
R(X»S)=(XK0L-X)**2+ALSQ*( < SKOL-S ) ♦♦a+NKL^^a ) 
IF(ABS(NKL/DERS).GT,»0001) 60 TO 10 
IF(RM(N) ,GT.1. ) GO TO 10 

IF( (ABS( <XKOL-X)/DERS) •LT.« 0001)* AND. (SKOL«LE«S))60T020 
10 IF(RM(N)«GT.l.) GO TO 12 

FUNKUsAL06<SQAL*(SKOL-S)<t>SQRT<K(X«S) ) ) 
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GO TO 30 

12 1F(UXK0L-X)/DERS).LT..0001) GO TO 20 
FUNKUsASXN( (SK0L">S ) OSQAL/ (XK0t*X ) ) 

GO TO 30 
20 FUNKUsO 
30 CONTINUE 
RETURN 
END 


SUBROUTINE CMATR(NlfNTOTWf NAD0WfNM0DEl*NM0DE2t NTEMP* 
lNSIZEl«NSIZE2«NCSIZEfNADDCf NPP«NAD0W2vNADDW3fNAD0WF) 
DIMENSION DATA(600 ) * RM (10 ) • FRA < 10 ) » NVECT ( 5 ) t VECT ( 5) 

REAL NKL»NLK«NKLNLK 
COMPLEX A(6130)fCK 
EQUIVALENCE (NOATA* DATA) 

COMMON/BLOCKA/NDATAC800 ) f A f NRP ♦ RM « FR A i DERS • OC » ALSO t SQAL 
COMMON/BLOCKB/XKoLf SKOLfLSYMtNPLANE,N,NKL 
COMMON/BLOCKC/NVECTiVECT 
MRITE(6«115) 

115 FORMATUX/ZIX^PART 3-OERIVE CMATRIX»//) 

NL1=NDATA(4) 

NL2=NDATA(8) 

IF (NMODE2.NE.O) GO TO SO 
DO 65 JsltNSIZEl 
65 A(NA00W2-t>J)sA(NA0DW4j) 

80 DO 90 J=ltNSIZE2 

90 A(NADDW3'»'J)::CMPLX(REAL(A(NAD0W2-t>J) ) t FRA ( N ) «AIMAG ( A 
l(NADOW2-i^J) ) ) 

DO 95 J=lfNSIZE2 
95 A(NADDWF-fJ)s:(0,f0«) 

CALL INVER(AfNRP) 

DO 200 MOOEs 1»NTEHP 
DO 200 LM=lfNRP 
DO 200 KM=liNRP 
J1=N1*(LM-1)+KM 
J2=NRP*( MODE-1 )+KM 
J3=NRP* ( MODE-1) +LM 

200 A(NADDWF4J3)=A(NA0DWF'*'J3)'*’A( Jl)>t'A(NADDW3>J2) 

IF(NPP .EQ,0) GO TO 217 
WRITE(6fll7) 

117 F0RMAT(2X*LISTING OF POTENTIAL DISCONTINUITIES*) 

CALL ALIST(A»NADDWF*NTEMPtNRP) 

217 CONTINUE 

DO 219 JsltNCSlZE 

219 A(NADDC4J)=(0. «0t > 

DO 600 J=l«NTEMP 
DO 600 I=lfNMODEl 
DO 600 L sl*NRP 
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DO 600 NTYPE=3fl5 

CALL TEST(NTYPE*L*NL.NADDJ) 

IF(NL.EQ*0) GO TO 600 
IF(NOATA(NADOJ) .NE.L) GO TO 600 
Ji=NADDW+NRP*( I"1 )+L 
J2=NADDWF+NRP»( J-D+L 
J3=NADDC+NTEMP*(I-l)-»-J 
DO 335 LIN=1.NL2 
NA0D2rN0ATA(5)+4»(LlN -1) 

IF (N0ATA(nA002) .EQ.L ) GO TO 355 
335 CONTINUE 
GO TO 600 

355 XL=:DATA(NA002 + 2) 

NSURF=N0ATA(NADD2+1) 

DO 445 LINl=lfNLl 
NADD1L=N0ATA(1)+9»(LIN1-1) 

NSURL=NDATA(NADD1U 
IF(NSURF.EQ.NSURL) GO TO 455 
445 CONTINUE 
GO TO 600 

455 NGE0MsNDATA(NADDlL+3) 

JTYPE=NTYPE-2 

GO T0(230«240«250«260«270«260 t290f 300f 3l0,3a04330t340« 


1350 ) • JTYPE 
230 CALLCCONIK 
GO TO 500 
240 CALLCC0N12( 
GO TO 500 
250 CALL CC0N5( 
60 TO 500 
260 CALL CC0N6( 
GO TO 500 
270 GO TO 250 
280 GO TO 260 
290 CALL C0LK9 
GO TO 600 
300 CALL CBLK9 


Jl, J2,J3«NGE0M*XLtL,NA0DJ«CK) 
Jl» J2f J3«NGE0Mf XL«L«NAOOJfCK) 
Jit J2« J3«NGE0n«XLfLfNA0DJfCK) 
J1«J2, J3«NGEOM,XL,L«NADDJtCK) 


GO TO 600 

310 CALLCCONIK Jl,J2«J3t NGEOM, XL tLtNADOJtCK) 
GO TO 500 

320 CALLCCON12( Jit J2f J3tNGEOntXLtL«NADOJtCK> 
GO TO 500 
330 CALL CBLK13 
GO TO 600 
340 CALL C0LK14 
GO TO 600 
350 CALL CBLK15 
GO TO 600 

500 IF<N0ATA(NA00JKNE.U 60 TO 600 
A( J3)sA( J3)-i>CK 
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600 continue 

WRITE(6*dOO) Nf Rn(N) •FRA(N) 

800 FORf1AT(lX//lX*CMATRlX FOR 

HX*MACH N0.a#fF9.5t* FRAs*tF9,5/) 
CALL ALIST(A,NAOOC,NKOOElfNTCMP) 

999 CONTINUE 
RETURN 
END 


SUBROUTINE CC0N5 ( J1 • J2 • J3 • NGEOM t XL t L « NADD j , CK ) 

DIMENSION DATA(800 ) #RM (10 ) t FRA ( 10 ) * VECT ( 5 ) • NVECT C 5 ) 

REAL NKL«NLK,NKLNLK 
COMPLEX A(6130)tCK 
EQuIVALENCE(NOATA»DATA) 

COMMON/BLOCKA/nOATA ( 800 ) « A»NRP«RM«FRA»DERSf 0C« ALSQf SQAL 
COM«ON/0LOCKB/XKoL,SKOL»LSYM,NPLANEfN»NKL 
COMMON/0LOCKC/NVECT*VECT 
IF(NDATA(NA0DJ) .NE,L) 60 TO 290 
XFORsVECT(l) 

XAFT=VECT(2) 

SINsVECT(3» 

SOUT=VECT(*f » 

ALPs(XAFT"XFOR)*(SOUT-SIN) 

CK=-( ALP*(REAL(A( Jl) } • ( 0 • , 1 • ) 4^FRA ( N S *A1 MAg ( A ( jl ) ) )«A( j2> ) 
1 ♦NGEOM 
290 continue 
RETURN 
ENTRY CC0N6 

IF(N0ATA(NA00J) .NE. L) GO TO 390 
XF0R=VECT(2> 

XAFT=VECT( 3) 

SINsVECT(4) 

SOUT=VECT(5) 

NBOXES=NVECT(25 

DXsXFOR-XL 

DELXSDX+ « XAFT»XF0R I ♦ . 5/NBOXES 

CK= ( (S0UT-’SINI*(DX^REAL(A( Jl) ) 

1 + AIMAG( A( jin )*A( J2S*CEXP(-(0.^1. )*FRA(NUDELXn*N6E0M 
390 CONTINUE 
RETURN 
END 
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subroutine: CCONIK J3«NGC0n,XL«L«NA00JtCK) 

OINtNSlON OATA(600 ) »RM ( 10 ) • FRA ( 10 ) t NVECT ( 5 ) • VCCT ( 5 ) 

REAL NKLtNLKtNKLNLK 
COMPLEX A(6130) *CK 
EQUIVALENCE (NOATA* DATA) 

COMMON/BLOCKA/NDATA(600 ) fA«NRP<RM«FRA«DERS«OCf ALSQiSQAL 
COMMON/BLOCKB/XKOLtSKOL»LSYM»NPLANEtNfNKL 
COHMON/BLOCKC/NVECTiVECT 
IF(N0ATA(NAD0J) .NE.L) GO TO 290 
XFORsVECTd) 

XAFT=VECT(2) 

SlNsVECT(3) 

S0UT=VECT(4) 

CKs(SOUT-SIN)*< (AIHAG(A(J1 ) )+(XFOR-XL)*REaL(A(J 1) ) )• 
lCEXP(-(0, •!. )*FRA(N)*<XFOR-XL> ) - < AIMAG ( A ( Ji ) ) + ( XAFT-XL ) 
2»REAL(A( Jl) ) )»CEXP(-(0. »1. )*FRA(N)*(XAFT*XL) ) ) *A I J2 ) *NGEOH 

290 CONTINUE 

RETURN 

ENTRY CC0N12 

IF(NOATA(NAODJ).NE* L) GO TO 390 
XF0RsVECT(2) 

SlNsVECT(4) 

S0UT=VECT(5) 

CKatSOUT-SIN)*J ( AIMAG( A ( Jl ) ) -^ ( XFOR-XL ) ♦REAL ( A ( Jl ) ) ) ♦ 
lNGEOM*CEXP(-(0.tl# ) ♦FRA ( N ) ♦ ( XFOR-XL ) ) )^A( J2) 

390 CONTINUE 

return 

END 


subroutine AlIST ( AtNAODf NR0WS«NC0LS) 
COMPLEX A(6130) 

WRITE (6»201) 

DO 200 JsltNROWS 
DO 200 KsliNCOLS 
ISaNADD-f C J-1 ) ♦NCOLS+K 
URITE(Gf202) JiKtA(IS) 

201 FORMAT (X^R0W^3X^C0L^2X^REAL^14X^ IMAGINARY* ) 

202_FORMAT(X2I4*2E1^.AJ' ^ _ 

200 CONTINUE 

return 

END 
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SUBROUTINE TEST (NTYPE«L»NL»NAObj» 

DIMENSION OATA(800 ) f Rli(lO) «FRA(X0) «NVECT(5) tVECT(S) 
COMPLEX A(6130) 

EQUIVALENCE(NDATAfDATA) " 

COMMON/BLOCKA/NOATACaOO ) «A«NRPtRMtFRA«DERS«OCt ALSOf SOAL 
common/blockc/nvect,vect 

NADDsNDAT A ( 4»NTYPE-3 ) 

NfIXsN0ATA(4*NTYPE-2) 

NFL=N0ATA(4*NTYPE«XJ _ 

NL=NDATA(4*NTYPE) 

IF(NL.EQ.O) go to 697 

DO 159 LlNslfNL 

NAODJs:NAOD+ ( LIN.l ) ♦ ( NFIX^NFL ) 

IF(NDATA(NA0DJ) .NE.L) go TO 159 

DO 157Jal,5 

NVECT(J)sO 
157 VECT(J)aO 

DO 156 JsliNPIX 
156 NVECT( J)=NDATA(NAODJ<fJ«l) 

NFLsNFL+NFlX-1 
DO 156 Jsl.NFL 
156 VECT( J)sOATA(NADDJ^J) 

GO TO 697 
159 CONTINUE 
897 RETURN 
END 
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SUBROUTINE INVER (A«N) 

COMPLEX PIVOT«A(6130) 

NlsN-M 
DO 1 Isl«N 
DO 1 JsltN 

1 A ( I*Nl-Nl+N+2-J)sAa»Nl*Nl+N+X»J) 

J=X 

1 = 0 

2 1=1+1 
JX=J+1 
ID=0 
NN=0 

3 IF(REAL(A(I*N1-NX+JX) ) ) 5tH,5 

4 IF(AIMAG(A(1«NX-NX+JX) n 5«6t5 

5 IF(IO) 6«X0«6 

6 IF(l-N) 8f7*8 

7 1 = 0 

8 I=I+X 
NN=NN4X 

IF(N-NN) 20f20f3 
XO I0=J 

DO XOl LA=X»N 
XOX A(LA+NX*NX+J) = (0»«0. ) 

A<I*NX«NX+J) = a»Of 0.0) 

PIVOT=A<I*NX-NX+JX) 

DO XX M=1#NX 

XX An*NX-NX+M)=A(I*NX-Nl+M)/PlVOT 
DO 17 LA=X#N 
IF(l-'LA)X2,X2tX5 

12 KK=LA+1 
IF(N-KK) X8fl3«X3 

13 PIVOT=A(KK*NX-NX+JX) 

DO 14 M=1*NX 

14 A(KK*NX-NX+M)=A<KK*NX-NX+M)-PIV0T*A(I*NX»NX+M) 
GO TO 17 

15 PIVOT=A(LA*NX-NX+JX) 

DO 16 Msl.NX 

16 A(LA*N1-N1+M)=A(LA+NX-NX+M)-PIV0T*A(I*NX-N1+M) 

17 CONTINUE 

18 IF(N-J) 20f20fX9 

19 J=JX 

GO TO 2 

20 RETURN 
END 
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List of Integral Formulas 


i. 


2 . 


3 . 


4 . 


5 . 


6 , 


/ 

/ 

/ 

/ 

/ 

/ 


dx 


(ax^ + c)V2 c^x^ + c 
dx - >4ix^ + c 


+ c 


cx 


X dx 


(ax2 + c)^/2 a>/ax2 


+ c 


, — — = log (x^^+-s/ax^ + c), a > 0 

^x2 + c va 

- 1 + a? 


dx 


(x^ + b2)>/<2 + - b2 x>^^ - b^ 


dx 


x-N^x^ + b 


= 7= '°3 { 


tan 

Vax^ + b -“N/b 




b > c 
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Subroutines For Steady Supersonic Calculations 

SUBROUTINE AC0N3 ( TIKL « AA * NKLNLK , AK ♦ XL » SL , L • K t NADOJ ) 
DIMENSION DATA( 800 ) « RM UO ) • FR A ( 10 ) t NVECT ( 5 ) * VECT ( 5 ) 

REAL NKLfNLK^NKLNLK 
COMPLEX A<6130) f AK,AA 
EQUI VALENCE (NDATA* DATA) 

COMMON/BLOCKA/NDATA (800 ) » A » NRP » RM ♦ FRA • DErs t OC t ALSO • SQAL 
COMMON/BLOCKB/XKOLtSKOL*LSYMtNPLANE«N«NKL 
COMMON/BLOCKC/NvECT« VECT 
R (XiS)=(XK0L-X)**2+ALSQ*( ( SKOL-S ) **2+NKL*#2 ) 

AA=0 

AKS::0 

IF(NDATA(NADDJ) .EQ.L) GO TO 100 
GO TO 894 
100 XFOR=VECT(l) 

XAFT=VECT(2> 

SIN=VECT(3) 

S0UT=VECT(4) 

NBOXES=l 
GO TO 267 
ENTRY ACON4 
AA=0 
AKS=0 

IF(NDATACNAODJ) •EQ.L) 60 TO 200 
GO TO 894 

200 NB0XES=NVECT(2) 

XF0R=VECT(2) 

XAFT=VECT(3) 

SIN=VECT(4) 

S0UT=VECT(5) 

267 DO 894 NB0X=1 • NBOXES 
500 IF(XL.GT.XKOL) GO TO 894 

IF(L,EQ.K.AND.SK0L,EQ.SL) go to 278 
IF{SL.EQ.SKOL.AND.XKOL.EQ.XAFT) GO TO 279 

277 AK=( 1./12. 56637) ♦ ( FUNS ( XAFT ♦ SOUT ) -FUNS ( XaFT t SIN ) -FUNS 
l(XFOR,SOUT)+FUNS(XFORfSIN) ) 

AK=(2.»A)5 + AKS)*LSYM*T1KL 
GO TO 884 

278 AKS=SQAL/(2.*( XAFT-XFOR) ) 
ak=aks*tikl*lsym 

GO TO 884 

279 AKS=-SQAL/(2.=«'(XAFT-XF0R) ) 

GO TO 277 

884 AA = AA-t>AK 
894 continue 
RETURN 
END 
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Cont’d Subroutines for Steady Supersonic Calculations 


FUNCTION FUNS(XfS) 

REAL NKLtNLKtNKLNLK 

DIMENSION OATA(800 ) * RM ( 10 ) » FRA { 10 ) » NVECT ( 5 ) , VECT ( 5 ) 
COMPLEX A( 6130) 

EQUlVALENCE(NDATAfDATA) 

C0MM0N/8L0CKA/NDATA( 800) t A t NRP • RM t FRA t DE rS * OC * ALSQ • SQAL 
COMMON/BLOCK0/XKOL,SKOL»LSYMiNPLANE«N*NKL 
R(X«S) = < XK0L-X)**2+ALSQ*( ( SKOL-S ) ♦*2+NKL**2 ) 
IF(NPlANE.N£.0 ) 60 TO 10 
RETURN 
10 CONTINUE 

IF(ABSUXK0L-X)/DERS),LT..0001 ) GO TO 20 
IF( ABS( (SK0L-S)/DERS) .LT, .0001) GO TO 20 
ARG=R( XtS) 

IF( -SQAL*ABS(S-SKOL) .LE. (X-XKOL) ) GO TO 2o 
IFiARG.GT.-. 000001 .AND. ARG.LT.O) ARG=fi 
FUNS=SQRT( AK6)/( ( XKOL-X ) ♦ ( SK0L-»S ) ) 

60 TO 30 
20 FUNS=0 
30 CONTINUE 

return 

END 
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Modifications In DMATR for Steady Supersonic Calculations: The method for 

breaking up a surface into rectangular elements in the case of steady 
supersonic calculation is slightly different from that used in subsonic case, 
and is shown In Fig. 36. The modifications required in subroutine DMATR 
are as fol lows: 

Subsonic 

DATA(K3)=SF0RCM)+(.75+(IQ-| ))* 

DELX(M) 

XF=SF0R(M)+DELX(M)^.25 
DATA(K5)=SAFT(M)-DELX(M)*.75 

The DMATR program listed in Appendix A is in the subsonic form. 
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